apply runic formatting

Author: briochemc

Project.toml

@@ -1,6 +1,6 @@
 name = "OceanTransportMatrixBuilder"
 uuid = "c2b4a04e-6049-4fc4-aa6a-5508a29a1e1c"
-version = "0.8.0"
+version = "0.8.1"
 authors = ["Benoit Pasquier <[email protected]> and contributors"]
 
 [deps]

src/RediGM.jl

@@ -1,6 +1,3 @@
-
-
-
 function localdensityslope(ρ, lon, lat, Z, I, gridtopology, dir)
     distances = verticalfacetriadgroupdistances(lon, lat, Z, I, gridtopology, dir)
     localρ = verticalfacetriadgroupvalues(ρ, I, gridtopology, dir)
@@ -38,9 +35,6 @@ function globalpotentialdensityslope(gsw_rho, so, ct, gridmetrics, indices, dir)
 end
 
 
-
-
-
 """
     bolus_GM_velocity(σ, gridmetrics; κGM = 600, maxslope = 0.01)
 
@@ -65,7 +59,7 @@ function bolus_GM_velocity(ρ, gridmetrics, indices; κGM = 600, maxslope = 0.01
     Sc = 0.004
     Sd = 0.001
     # that should not work since Sᵢ and Sⱼ are not colocated
-    taper = @. 0.5 * (1 + tanh((Sc - sqrt(Sᵢ ^ 2 + Sⱼ ^ 2)) / Sd))
+    taper = @. 0.5 * (1 + tanh((Sc - sqrt(Sᵢ^2 + Sⱼ^2)) / Sd))
     Sᵢ = taper .* Sᵢ
     Sⱼ = taper .* Sⱼ
 

src/classicderivatives.jl

@@ -1,4 +1,3 @@
-
 ####################################
 # Forward and backward derivatives #
 ####################################
@@ -78,10 +77,6 @@ function ∂ₖ₋(χ, gridmetrics)
 end
 
 
-
-
-
-
 # TODO: Check that I am properly dealing with unevenly spaced grids
 # The commented code below could help with that, as it provides the weights
 # for a finite difference stencil for any arrangement of grid points
@@ -100,5 +95,3 @@ end
 #     A = @. (x' - x₀)^ℓ / factorial(ℓ)
 #     return A \ (ℓ .== m) # vector of weights w
 # end
-
-

src/dyads.jl

@@ -1,5 +1,3 @@
-
-
 #########
 # DYADS #
 #########
@@ -78,4 +76,3 @@ function globalverticaldyadderivative(χ, gridmetrics, indices)
     end
     return ∂χ
 end
-

src/extratools.jl

@@ -35,7 +35,7 @@ where the lumping should occur.
 Outside of this region, no lumping.
 Default is `f=Returns(true)`, i.e. lump everywhere.
 """
-function lump_and_spray(wet3D, vol, T, mask=trues(size(wet3D)); di=2, dj=2, dk=1)
+function lump_and_spray(wet3D, vol, T, mask = trues(size(wet3D)); di = 2, dj = 2, dk = 1)
 
     # extend the grid to avoid lumping cells outside of bounds
     nxyz = size(wet3D)
@@ -53,15 +53,15 @@ function lump_and_spray(wet3D, vol, T, mask=trues(size(wet3D)); di=2, dj=2, dk=1
 
     # loop once over all grid cells and assign lumped indices
     c = 2 # index / counter (we start at 2 because 1 is reserved for dry cells)
-    neighbours = CartesianIndices((0:di-1, 0:dj-1, 0:dk-1))
+    neighbours = CartesianIndices((0:(di - 1), 0:(dj - 1), 0:(dk - 1)))
     for 𝑖 in eachindex(C)
         C𝑖 = C[𝑖]
         # skip if index already assigned and in mask
         #(if assigned but not in mask, must reassign)
         LUMPidx[C𝑖] > 0 && mask[C𝑖] && continue
         if mask[C𝑖] # if in region, assign all neighbours also in region
             # list of neighbours
-            L𝑖 = L[C𝑖 .+ neighbours][:]
+            L𝑖 = vec(L[C𝑖 .+ neighbours])
             # First, assign dry index to dry cells
             localwet = wet3Dext[L𝑖]
             dryidx = L𝑖[.!localwet]
@@ -105,23 +105,22 @@ function lump_and_spray(wet3D, vol, T, mask=trues(size(wet3D)); di=2, dj=2, dk=1
     nwet = size(LUMP, 2)
     nwet_c = size(LUMP, 1)
     @info """LUMP and SPRAY:
-    Matrix size reduction: $(round(100(1 - nwet_c/nwet)))% ($nwet -> $nwet_c)
+    Matrix size reduction: $(round(100(1 - nwet_c / nwet)))% ($nwet -> $nwet_c)
     """
 
     return LUMP, SPRAY, vol_c
 end
 
 
-
 """
 as2D(x, wet3D)
 
 `x` must be a vector with `length(x) == sum(wet3D[:,:,1])`
 """
 function as2D(x, wet3D)
     x2D = fill(NaN, size(wet3D)[1:2])
-    @views x2D[wet3D[:,:,1]] .= x
-    x2D
+    @views x2D[wet3D[:, :, 1]] .= x
+    return x2D
 end
 
 """
@@ -132,26 +131,5 @@ as3D(x, wet3D)
 function as3D(x, wet3D)
     x3D = fill(NaN, size(wet3D))
     @views x3D[wet3D] .= x
-    x3D
+    return x3D
 end
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-

src/gridcellgeometry.jl

@@ -1,5 +1,3 @@
-
-
 abstract type ArakawaGridCell end
 
 struct AGridCell <: ArakawaGridCell
@@ -70,8 +68,8 @@ function getarakawagrid(u_lon, u_lat, v_lon, v_lat, gridmetrics)
 
     cell = (; C, SW, SE, NE, NW, S, N, W, E)
 
-    u_distances = (; (k => haversine(P, u_point) for (k,P) in pairs(cell))...)
-    v_distances = (; (k => haversine(P, v_point) for (k,P) in pairs(cell))...)
+    u_distances = (; (k => haversine(P, u_point) for (k, P) in pairs(cell))...)
+    v_distances = (; (k => haversine(P, v_point) for (k, P) in pairs(cell))...)
 
     u_distance, u_pos = findmin(u_distances)
     v_distance, v_pos = findmin(v_distances)
@@ -126,8 +124,8 @@ function interpolateontodefaultCgrid(u, u_lon, u_lat, v, v_lon, v_lat, gridmetri
     # so that's what we do here
     u2 = replace(u, _FillValue => 0.0)
     v2 = replace(v, _FillValue => 0.0)
-    u2 = 0.5(u2 + [fill(0.0, nx, 1, nz);; u2[:, 1:end-1, :]])
-    v2 = 0.5(v2 + [fill(0.0, 1, ny, nz); v2[1:end-1, :, :]])
+    u2 = 0.5(u2 + [fill(0.0, nx, 1, nz);; u2[:, 1:(end - 1), :]])
+    v2 = 0.5(v2 + [fill(0.0, 1, ny, nz); v2[1:(end - 1), :, :]])
     SE_points = [(lon, lat) for (lon, lat) in zip(lon_vertices[2, :, :], lat_vertices[2, :, :])]
     NE_points = [(lon, lat) for (lon, lat) in zip(lon_vertices[3, :, :], lat_vertices[3, :, :])]
     NW_points = [(lon, lat) for (lon, lat) in zip(lon_vertices[4, :, :], lat_vertices[4, :, :])]
@@ -142,9 +140,6 @@ function interpolateontodefaultCgrid(u, u_lon, u_lat, v, v_lon, v_lat, gridmetri
 end
 
 
-
-
-
 """
     vertexpermutation(lon_vertices, lat_vertices)
 
@@ -167,8 +162,8 @@ function vertexpermutation(lon_vertices, lat_vertices)
     i = j = 1
     # Turn the vertices into points
     points = collect(zip(lon_vertices[:, i, j], lat_vertices[:, i, j]))
-    points_east = collect(zip(lon_vertices[:, i+1, j], lat_vertices[:, i+1, j]))
-    points_north = collect(zip(lon_vertices[:, i, j+1], lat_vertices[:, i, j+1]))
+    points_east = collect(zip(lon_vertices[:, i + 1, j], lat_vertices[:, i + 1, j]))
+    points_north = collect(zip(lon_vertices[:, i, j + 1], lat_vertices[:, i, j + 1]))
     # Find the common points
     common_east = Set(points) ∩ Set(points_east)
     common_noth = Set(points) ∩ Set(points_north)
@@ -183,8 +178,6 @@ function vertexpermutation(lon_vertices, lat_vertices)
 end
 
 
-
-
 # distances
 function horizontaldistance(lon, lat, I, J)
     I = horizontalindex(I)
@@ -207,42 +200,40 @@ verticalindex(I::CartesianIndex{1}) = I
 verticalindex(I::CartesianIndex{3}) = CartesianIndex(I.I[3])
 
 
-
-
 # The default orientation is the following:
 #
 #     4 ────┐ 3
 #           │
 #     1 ────┘ 2
 #
 function vertexindices(dir)
-    dir == :south ? (1, 2) :
-    dir == :east ? (2, 3) :
-    dir == :north ? (3, 4) :
-    dir == :west ? (1, 4) :
-    error()
+    return dir == :south ? (1, 2) :
+        dir == :east ? (2, 3) :
+        dir == :north ? (3, 4) :
+        dir == :west ? (1, 4) :
+        error()
 end
-vertexpoint(vlon, vlat, i, j, vertexidx) = (vlon[vertexidx,i,j], vlat[vertexidx,i,j])
+vertexpoint(vlon, vlat, i, j, vertexidx) = (vlon[vertexidx, i, j], vlat[vertexidx, i, j])
 function verticalfacewidth(vlon, vlat, i, j, dir)
     a, b = vertexindices(dir)
     A = vertexpoint(vlon, vlat, i, j, a)
     B = vertexpoint(vlon, vlat, i, j, b)
-    haversine(A, B)
+    return haversine(A, B)
 end
 verticalfacewidth(edge_length_2D, i, j, dir) = edge_length_2D[dir][i, j]
 
 function verticalfacearea(vlon, vlat, lev_bnds_or_thkcello, i, j, k, dir)
     height = cellthickness(lev_bnds_or_thkcello, i, j, k)
     width = verticalfacewidth(vlon, vlat, i, j, dir)
-    height * width
+    return height * width
 end
 function verticalfacearea(edge_length_2D, lev_bnds_or_thkcello, i, j, k, dir)
     height = cellthickness(lev_bnds_or_thkcello, i, j, k)
     width = verticalfacewidth(edge_length_2D, i, j, dir)
-    height * width
+    return height * width
 end
 
-cellthickness(lev_bnds::Matrix, i, j, k) = abs(lev_bnds[2,k] - lev_bnds[1,k])
+cellthickness(lev_bnds::Matrix, i, j, k) = abs(lev_bnds[2, k] - lev_bnds[1, k])
 cellthickness(thkcello, i, j, k) = thkcello[i, j, k]
 
 
@@ -252,11 +243,11 @@ function centroid2edgedistance(lon, lat, vlon, vlat, i, j, dir)
     A = vertexpoint(vlon, vlat, i, j, a)
     B = vertexpoint(vlon, vlat, i, j, b)
     M = midpointonsphere(A, B)
-    haversine(C, M)
+    return haversine(C, M)
 end
 
 function midpointonsphere(A, B)
-    if abs(A[1] - B[1]) < 180
+    return if abs(A[1] - B[1]) < 180
         (A .+ B) ./ 2
     else # if the edge crosses the longitudinal edge of the map
         (A .+ B) ./ 2 .+ (180, 0)
@@ -271,8 +262,6 @@ function horizontalcentroiddistance(lon, lat, iA, jA, iB, jB)
 end
 
 
-
-
 function makegridmetrics(; areacello, volcello, lon, lat, lev, lon_vertices, lat_vertices)
 
     # Differences in output "missing" data is commonplace,
@@ -282,19 +271,19 @@ function makegridmetrics(; areacello, volcello, lon, lat, lev, lon_vertices, lat
     haskey(volcello.properties, "_FillValue") && push!(toreplace, volcello.properties["_FillValue"])
     replacelist = (x => NaN for x in toreplace)
 
-	# volume (3D)
+    # volume (3D)
     v3D = volcello |> Array{Union{Missing, Float64}}
-	v3D = replace(v3D, replacelist...)
+    v3D = replace(v3D, replacelist...)
 
     # area (2D)
     area2D = areacello |> Array{Union{Missing, Float64}}
     area2D = replace(area2D, replacelist...)
 
-	# depth and cell height (3D)
-	thkcello = v3D ./ area2D
+    # depth and cell height (3D)
+    thkcello = v3D ./ area2D
     ZBOT3D = cumsum(thkcello, dims = 3)
     Z3D = ZBOT3D - 0.5 * thkcello
-	zt = lev |> Array
+    zt = lev |> Array
 
     lat = lat |> Array
     lon = lon |> Array
@@ -305,19 +294,18 @@ function makegridmetrics(; areacello, volcello, lon, lat, lev, lon_vertices, lat
 
     # sort the vertices to match the default orientation
     vertexidx = vertexpermutation(lon_vertices, lat_vertices)
-    lon_vertices = lon_vertices[vertexidx,:,:]
-    lat_vertices = lat_vertices[vertexidx,:,:]
+    lon_vertices = lon_vertices[vertexidx, :, :]
+    lat_vertices = lat_vertices[vertexidx, :, :]
 
     C = CartesianIndices(size(lon))
 
     gridtopology = getgridtopology(lon_vertices, lat_vertices, zt)
 
     dirs = (:south, :east, :north, :west)
     𝑗s = (j₋₁, i₊₁, j₊₁, i₋₁)
-    edge_length_2D = Dict(d=>[verticalfacewidth(lon_vertices, lat_vertices, 𝑖.I[1], 𝑖.I[2], d) for 𝑖 in C] for d in dirs)
-    distance_to_edge_2D = Dict(d=>[centroid2edgedistance(lon, lat, lon_vertices, lat_vertices, 𝑖.I[1], 𝑖.I[2], d) for 𝑖 in C] for d in dirs)
-    distance_to_neighbour_2D = Dict(d=>[horizontaldistance(lon, lat, 𝑖, 𝑗(𝑖, gridtopology)) for 𝑖 in C] for (d,𝑗) in zip(dirs,𝑗s))
+    edge_length_2D = Dict(d => [verticalfacewidth(lon_vertices, lat_vertices, 𝑖.I[1], 𝑖.I[2], d) for 𝑖 in C] for d in dirs)
+    distance_to_edge_2D = Dict(d => [centroid2edgedistance(lon, lat, lon_vertices, lat_vertices, 𝑖.I[1], 𝑖.I[2], d) for 𝑖 in C] for d in dirs)
+    distance_to_neighbour_2D = Dict(d => [horizontaldistance(lon, lat, 𝑖, 𝑗(𝑖, gridtopology)) for 𝑖 in C] for (d, 𝑗) in zip(dirs, 𝑗s))
 
-	return (; area2D, v3D, thkcello, lon_vertices, lat_vertices, lon, lat, Z3D, zt, edge_length_2D, distance_to_edge_2D, distance_to_neighbour_2D, gridtopology)
+    return (; area2D, v3D, thkcello, lon_vertices, lat_vertices, lon, lat, Z3D, zt, edge_length_2D, distance_to_edge_2D, distance_to_neighbour_2D, gridtopology)
 end
-

src/gridtopology.jl

@@ -25,16 +25,16 @@ function getgridtopology(lon_vertices, lat_vertices, lev)
     ny = size(lon_vertices, 3)
     nz = length(lev)
     # North pole vertices
-    NPlon = @view lon_vertices[3:4,:,end]
-    NPlat = @view lat_vertices[3:4,:,end]
+    NPlon = @view lon_vertices[3:4, :, end]
+    NPlat = @view lat_vertices[3:4, :, end]
     # If all the latitudes of the "northmost" vertices are 90, then it's a "regular" grid with 2 poles
     if all(NPlat .== 90)
-        return BipolarGridTopology(nx,ny,nz)
-    # Otherwise check if the north pole is split in two
+        return BipolarGridTopology(nx, ny, nz)
+        # Otherwise check if the north pole is split in two
     elseif isapprox(NPlon, rot180(NPlon)) && isapprox(NPlat, rot180(NPlat))
-        return TripolarGridTopology(nx,ny,nz)
+        return TripolarGridTopology(nx, ny, nz)
     else
-        return UnknownGridTopology(nx,ny,nz)
+        return UnknownGridTopology(nx, ny, nz)
     end
 end
 
@@ -55,20 +55,19 @@ k₋₁(C, ::AbstractGridTopology) = C.I[3] > 1 ? C + CartesianIndex(0, 0, -1) :
 getindexornan(χ, I) = isnothing(I) ? NaN : χ[I] # TODO check that this is going to work
 
 
-
-function ishift(C, g::AbstractGridTopology, n=0)
+function ishift(C, g::AbstractGridTopology, n = 0)
     i, j, k = C.I
-    CartesianIndex(mod1(i + n, g.nx), j, k)
+    return CartesianIndex(mod1(i + n, g.nx), j, k)
 end
-function jshift(C, g::AbstractGridTopology, n=0)
+function jshift(C, g::AbstractGridTopology, n = 0)
     i, j, k = C.I
     j2 = j + n
-    (1 ≤ j2 ≤ g.ny) ? CartesianIndex(i, j2, k) : nothing
+    return (1 ≤ j2 ≤ g.ny) ? CartesianIndex(i, j2, k) : nothing
 end
-function kshift(C, g::AbstractGridTopology, n=0)
+function kshift(C, g::AbstractGridTopology, n = 0)
     i, j, k = C.I
     k2 = k + n
-    (1 ≤ k2 ≤ g.nz) ? CartesianIndex(i, j, k2) : nothing
+    return (1 ≤ k2 ≤ g.nz) ? CartesianIndex(i, j, k2) : nothing
 end
 
 # Special behavior for tripolar grids where the seam connects the top and "folds" it around
@@ -81,7 +80,7 @@ end
 j₊₁(C, g::TripolarGridTopology) = C.I[2] < g.ny ? C + CartesianIndex(0, 1, 0) : CartesianIndex(g.nx - C.I[1] + 1, g.ny, C.I[3])
 j₊₁(C::CartesianIndex{2}, g::TripolarGridTopology) = C.I[2] < g.ny ? C + CartesianIndex(0, 1) : CartesianIndex(g.nx - C.I[1] + 1, g.ny)
 
-function jshift(C, g::TripolarGridTopology, n=0)
+function jshift(C, g::TripolarGridTopology, n = 0)
     i, j, k = C.I
     j2 = j + n
     if !(1 ≤ j2)