From 1a581e6f9ce7475549056500f4e1c56a8335ca77 Mon Sep 17 00:00:00 2001 From: Benoit Pasquier Date: Thu, 10 Oct 2024 20:52:05 +1100 Subject: [PATCH] Add new topology --- Project.toml | 2 +- src/OceanTransportMatrixBuilder.jl | 3 +- src/{gridtypes.jl => grids.jl} | 35 ++++- src/matrixbuilding.jl | 230 +++++++++++++++-------------- src/preprocessing.jl | 32 +--- src/topology.jl | 56 +++++++ test/localbuild.jl | 49 +++++- 7 files changed, 257 insertions(+), 150 deletions(-) rename src/{gridtypes.jl => grids.jl} (76%) create mode 100644 src/topology.jl diff --git a/Project.toml b/Project.toml index 079f599..278bca2 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "OceanTransportMatrixBuilder" uuid = "c2b4a04e-6049-4fc4-aa6a-5508a29a1e1c" authors = ["Benoit Pasquier and contributors"] -version = "0.2.2" +version = "0.2.3" [deps] Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7" diff --git a/src/OceanTransportMatrixBuilder.jl b/src/OceanTransportMatrixBuilder.jl index 766a635..e6d006b 100644 --- a/src/OceanTransportMatrixBuilder.jl +++ b/src/OceanTransportMatrixBuilder.jl @@ -9,7 +9,8 @@ using Distances: haversine include("preprocessing.jl") include("matrixbuilding.jl") -include("gridtypes.jl") +include("grids.jl") +include("topology.jl") include("extratools.jl") # <- I think this should be in a separate "base" repo export velocity2fluxes diff --git a/src/gridtypes.jl b/src/grids.jl similarity index 76% rename from src/gridtypes.jl rename to src/grids.jl index 4a8c289..2adcd47 100644 --- a/src/gridtypes.jl +++ b/src/grids.jl @@ -1,3 +1,35 @@ +# The default orientation is the following: +# +# 4 ────┐ 3 +# │ +# 1 ────┘ 2 +# +# Given lon_vertices and lat_vertices, find the permutation +# that sorts the vertices in that order. +function vertexpermutation(lon_vertices, lat_vertices) + # Make sure the vertices are in the right shape (4, nx, ny) + @assert size(lon_vertices, 1) == size(lat_vertices, 1) == 4 + # Take the first grid cell + 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])) + # Find the common points + common_east = Set(points) ∩ Set(points_east) + common_noth = Set(points) ∩ Set(points_north) + # Find the indices of the common points + idx_east = findall(in(common_east), points) + idx_north = findall(in(common_noth), points) + idx3 = only(idx_east ∩ idx_north) # common to all 3 cells + idx2 = only(setdiff(idx_east, idx3)) # common to (i,j) and (i+1,j) only + idx4 = only(setdiff(idx_north, idx3)) # common to (i,j) and (i,j+1) only + idx1 = only(setdiff(1:4, idx2, idx3, idx4)) # only in (i,j) + return [idx1, idx2, idx3, idx4] +end + + + """ arakawa, uo_pos, vo_pos, uo_relerr, vo_relerr = gridtype(uo_lon, uo_lat, vo_lon, vo_lat, modelgrid) @@ -114,4 +146,5 @@ function interpolateontodefaultCgrid(uo, uo_lon, uo_lat, vo, vo_lon, vo_lat, mod error("Interpolation not implemented for this grid type") end -end \ No newline at end of file +end + diff --git a/src/matrixbuilding.jl b/src/matrixbuilding.jl index fdceeef..2e8dbfe 100644 --- a/src/matrixbuilding.jl +++ b/src/matrixbuilding.jl @@ -131,7 +131,7 @@ end # There are three ways to index: # - Cartesian indices (i,j,k) -# - Linear index Li +# - Linear index L𝑖 # - Wet linear index 𝑖 (that's the one we want to record for the matrix) # So to fill T[𝑖,𝑗] -ϕ[𝑖→𝑗] / m[𝑖], I need be able to convert, in sequence: # 𝑖 -> (i,j,k) -> neihghbour (i′,j′,k′) -> 𝑗 @@ -146,7 +146,7 @@ Return the sparse (i, j, v) for the upwind advection operator Tadv. function upwind_advection_operator_sparse_entries(; ϕ, modelgrid, indices, ρ) # Unpack model grid - (; v3D,) = modelgrid + (; v3D, gridtype) = modelgrid # Unpack indices (; wet3D, Lwet, Lwet3D, C) = indices @@ -156,62 +156,56 @@ function upwind_advection_operator_sparse_entries(; ϕ, modelgrid, indices, ρ) nx, ny, _ = nxyz @time for 𝑖 in eachindex(Lwet) - Li = Lwet[𝑖] - i, j, k = C[Li].I - m𝑖 = v3D[i,j,k] * ρ + L𝑖 = Lwet[𝑖] + C𝑖 = C[L𝑖] + i, j, k = C𝑖.I + m𝑖 = v3D[C𝑖] * ρ # From West - ϕwest = ϕ.west[i,j,k] + ϕwest = ϕ.west[C𝑖] if ϕwest > 0 - i′ = mod1(i - 1, nx) - 𝑗 = Lwet3D[i′,j,k] - ismissing(𝑗) && @show(i, j, k, i′) - m𝑗 = v3D[i′,j,k] * ρ + C𝑗 = i₋₁(C𝑖, gridtype) + 𝑗 = Lwet3D[C𝑗] + m𝑗 = v3D[C𝑗] * ρ pushTadvectionvalues!(𝑖s, 𝑗s, Tvals, 𝑖, 𝑗, ϕwest, m𝑖, m𝑗) end # From East - ϕeast = ϕ.east[i,j,k] + ϕeast = ϕ.east[C𝑖] if ϕeast < 0 - i′ = mod1(i + 1, nx) - 𝑗 = Lwet3D[i′,j,k] - m𝑗 = v3D[i′,j,k] * ρ + C𝑗 = i₊₁(C𝑖, gridtype) + 𝑗 = Lwet3D[C𝑗] + m𝑗 = v3D[C𝑗] * ρ pushTadvectionvalues!(𝑖s, 𝑗s, Tvals, 𝑖, 𝑗, -ϕeast, m𝑖, m𝑗) end # From South - ϕsouth = ϕ.south[i,j,k] + ϕsouth = ϕ.south[C𝑖] if ϕsouth > 0 - j′ = j - 1 - 𝑗 = Lwet3D[i,j′,k] - m𝑗 = v3D[i,j′,k] * ρ + C𝑗 = j₋₁(C𝑖, gridtype) + 𝑗 = Lwet3D[C𝑗] + m𝑗 = v3D[C𝑗] * ρ pushTadvectionvalues!(𝑖s, 𝑗s, Tvals, 𝑖, 𝑗, ϕsouth, m𝑖, m𝑗) end # From North (Special case with north bipole) - ϕnorth = ϕ.north[i,j,k] + ϕnorth = ϕ.north[C𝑖] if ϕnorth < 0 - if j == ny - j′ = j - i′ = nx - i + 1 - else - j′ = j + 1 - i′ = i - end - 𝑗 = Lwet3D[i′,j′,k] - m𝑗 = v3D[i′,j′,k] * ρ + C𝑗 = j₊₁(C𝑖, gridtype) + 𝑗 = Lwet3D[C𝑗] + m𝑗 = v3D[C𝑗] * ρ pushTadvectionvalues!(𝑖s, 𝑗s, Tvals, 𝑖, 𝑗, -ϕnorth, m𝑖, m𝑗) end # From Bottom - ϕbottom = ϕ.bottom[i,j,k] + ϕbottom = ϕ.bottom[C𝑖] if ϕbottom > 0 - k′ = k + 1 - 𝑗 = Lwet3D[i,j,k′] - m𝑗 = v3D[i,j,k′] * ρ + C𝑗 = k₊₁(C𝑖, gridtype) + 𝑗 = Lwet3D[C𝑗] + m𝑗 = v3D[C𝑗] * ρ pushTadvectionvalues!(𝑖s, 𝑗s, Tvals, 𝑖, 𝑗, ϕbottom, m𝑖, m𝑗) end # From Top - ϕtop = ϕ.top[i,j,k] + ϕtop = ϕ.top[C𝑖] if ϕtop < 0 && k > 1 # Evaporation/precipitation -> no change to χ - k′ = k - 1 - 𝑗 = Lwet3D[i,j,k′] - m𝑗 = v3D[i,j,k′] * ρ + C𝑗 = k₋₁(C𝑖, gridtype) + 𝑗 = Lwet3D[C𝑗] + m𝑗 = v3D[C𝑗] * ρ pushTadvectionvalues!(𝑖s, 𝑗s, Tvals, 𝑖, 𝑗, -ϕtop, m𝑖, m𝑗) end end @@ -276,7 +270,7 @@ Return the sparse (i, j, v) for the horizontal diffusion operator TκH. function horizontal_diffusion_operator_sparse_entries(; modelgrid, indices, κH, ΩH) # Unpack model grid - (; v3D, edge_length_2D, distance_to_edge_2D, DZT3d) = modelgrid + (; v3D, edge_length_2D, distance_to_edge_2D, DZT3d, gridtype) = modelgrid # Unpack indices (; wet3D, Lwet, Lwet3D, C) = indices @@ -286,65 +280,78 @@ function horizontal_diffusion_operator_sparse_entries(; modelgrid, indices, κH, @time for 𝑖 in eachindex(Lwet) ΩH[𝑖] || continue # only continue if inside ΩH - Li = Lwet[𝑖] - i, j, k = C[Li].I - V = v3D[i,j,k] + L𝑖 = Lwet[𝑖] + C𝑖 = C[L𝑖] + i, j, k = C𝑖.I + V = v3D[C𝑖] # From West - iW, jW = mod1(i - 1, nx), j - 𝑗W = Lwet3D[iW, jW, k] - # (𝑖 == 𝑗W) && @show(i, j, iW, jW) - if !ismissing(𝑗W) && ΩH[𝑗W] - # I take the minimum area from both dirs (through which mixing goes through) - aij = verticalfacearea(edge_length_2D, DZT3d, i, j, k, :west) - aji = verticalfacearea(edge_length_2D, DZT3d, iW, jW, k, :east) - a = min(aij, aji) - # I take the mean distance from both dirs - dij = horizontalcentroiddistance(distance_to_edge_2D, i, j, iW, jW, :west) - dji = horizontalcentroiddistance(distance_to_edge_2D, iW, jW, i, j, :east) - d = (dij + dji) / 2 - pushTmixingvalues!(𝑖s, 𝑗s, Tvals, 𝑖, 𝑗W, κH, a, d, V) + C𝑗W = i₋₁(C𝑖, gridtype) + if !isnothing(C𝑗W) + 𝑗W = Lwet3D[C𝑗W] + if !ismissing(𝑗W) && ΩH[𝑗W] + iW, jW, _ = C𝑗W.I + # (𝑖 == 𝑗W) && @show(i, j, iW, jW) + # I take the minimum area from both dirs (through which mixing goes through) + aij = verticalfacearea(edge_length_2D, DZT3d, i, j, k, :west) + aji = verticalfacearea(edge_length_2D, DZT3d, iW, jW, k, :east) + a = min(aij, aji) + # I take the mean distance from both dirs + dij = horizontalcentroiddistance(distance_to_edge_2D, i, j, iW, jW, :west) + dji = horizontalcentroiddistance(distance_to_edge_2D, iW, jW, i, j, :east) + d = (dij + dji) / 2 + pushTmixingvalues!(𝑖s, 𝑗s, Tvals, 𝑖, 𝑗W, κH, a, d, V) + end end # From East - iE, jE = mod1(i + 1, nx), j - 𝑗E = Lwet3D[iE, jE, k] - # (𝑖 == 𝑗E) && @show(i, j, iE, jE) - if !ismissing(𝑗E) && ΩH[𝑗E] - aij = verticalfacearea(edge_length_2D, DZT3d, i, j, k, :east) - aji = verticalfacearea(edge_length_2D, DZT3d, iE, jE, k, :west) - a = min(aij, aji) - dij = horizontalcentroiddistance(distance_to_edge_2D, i, j, iE, jE, :east) - dji = horizontalcentroiddistance(distance_to_edge_2D, iE, jE, i, j, :west) - d = (dij + dji) / 2 - pushTmixingvalues!(𝑖s, 𝑗s, Tvals, 𝑖, 𝑗E, κH, a, d, V) + C𝑗E = i₊₁(C𝑖, gridtype) + if !isnothing(C𝑗E) + 𝑗E = Lwet3D[C𝑗E] + if !ismissing(𝑗E) && ΩH[𝑗E] + iE, jE, _ = C𝑗E.I + # (𝑖 == 𝑗E) && @show(i, j, iE, jE) + aij = verticalfacearea(edge_length_2D, DZT3d, i, j, k, :east) + aji = verticalfacearea(edge_length_2D, DZT3d, iE, jE, k, :west) + a = min(aij, aji) + dij = horizontalcentroiddistance(distance_to_edge_2D, i, j, iE, jE, :east) + dji = horizontalcentroiddistance(distance_to_edge_2D, iE, jE, i, j, :west) + d = (dij + dji) / 2 + pushTmixingvalues!(𝑖s, 𝑗s, Tvals, 𝑖, 𝑗E, κH, a, d, V) + end end # From South - if j > 1 - iS, jS = i, j - 1 - 𝑗S = Lwet3D[iS, jS, k] - # (𝑖 == 𝑗S) && @show(i, j, iS, jS) - if !ismissing(𝑗S) && ΩH[𝑗S] - aij = verticalfacearea(edge_length_2D, DZT3d, i, j, k, :south) - aji = verticalfacearea(edge_length_2D, DZT3d, iS, jS, k, :north) - a = min(aij, aji) - dij = horizontalcentroiddistance(distance_to_edge_2D, i, j, iS, jS, :south) - dji = horizontalcentroiddistance(distance_to_edge_2D, iS, jS, i, j, :north) - d = (dij + dji) / 2 - pushTmixingvalues!(𝑖s, 𝑗s, Tvals, 𝑖, 𝑗S, κH, a, d, V) - end - end + C𝑗S = j₋₁(C𝑖, gridtype) + if !isnothing(C𝑗S) + 𝑗S = Lwet3D[C𝑗S] + if !ismissing(𝑗S) && ΩH[𝑗S] + iS, jS, _ = C𝑗S.I + # (𝑖 == 𝑗S) && @show(i, j, iS, jS) + aij = verticalfacearea(edge_length_2D, DZT3d, i, j, k, :south) + aji = verticalfacearea(edge_length_2D, DZT3d, iS, jS, k, :north) + a = min(aij, aji) + dij = horizontalcentroiddistance(distance_to_edge_2D, i, j, iS, jS, :south) + dji = horizontalcentroiddistance(distance_to_edge_2D, iS, jS, i, j, :north) + d = (dij + dji) / 2 + pushTmixingvalues!(𝑖s, 𝑗s, Tvals, 𝑖, 𝑗S, κH, a, d, V) + end + end # From North - # Note that the opposite direction (oppdir) is still north at j == ny - (iN, jN, oppdir) = (j == ny) ? (nx - i + 1, j, :north) : (i, j + 1, :south) - 𝑗N = Lwet3D[iN, jN, k] - # (𝑖 == 𝑗N) && @show(i, j, iN, jN) - if !ismissing(𝑗N) && ΩH[𝑗N] - aij = verticalfacearea(edge_length_2D, DZT3d, i, j, k, :north) - aji = verticalfacearea(edge_length_2D, DZT3d, iN, jN, k, oppdir) - a = min(aij, aji) - dij = horizontalcentroiddistance(distance_to_edge_2D, i, j, iN, jN, :north) - dji = horizontalcentroiddistance(distance_to_edge_2D, iN, jN, i, j, oppdir) - d = (dij + dji) / 2 - pushTmixingvalues!(𝑖s, 𝑗s, Tvals, 𝑖, 𝑗N, κH, a, d, V) + C𝑗N = j₊₁(C𝑖, gridtype) + if !isnothing(C𝑗N) + 𝑗N = Lwet3D[C𝑗N] + if !ismissing(𝑗N) && ΩH[𝑗N] + # (𝑖 == 𝑗N) && @show(i, j, iN, jN) + iN, jN, _ = C𝑗N.I + # Note that the opposite direction (oppdir) is still north at j == ny + # TODO: implement this into a topology.jl function + oppdir = (j == ny) ? :north : :south + aij = verticalfacearea(edge_length_2D, DZT3d, i, j, k, :north) + aji = verticalfacearea(edge_length_2D, DZT3d, iN, jN, k, oppdir) + a = min(aij, aji) + dij = horizontalcentroiddistance(distance_to_edge_2D, i, j, iN, jN, :north) + dji = horizontalcentroiddistance(distance_to_edge_2D, iN, jN, i, j, oppdir) + d = (dij + dji) / 2 + pushTmixingvalues!(𝑖s, 𝑗s, Tvals, 𝑖, 𝑗N, κH, a, d, V) + end end end @@ -378,7 +385,7 @@ end function vertical_diffusion_operator_sparse_entries(; modelgrid, indices, κV, Ω) # Unpack model grid - (; v3D, area2D, zt) = modelgrid + (; v3D, area2D, zt, gridtype) = modelgrid # Unpack indices (; wet3D, Lwet, Lwet3D, C) = indices @@ -388,28 +395,31 @@ function vertical_diffusion_operator_sparse_entries(; modelgrid, indices, κV, @time for 𝑖 in eachindex(Lwet) Ω[𝑖] || continue # only continue if inside Ω - Li = Lwet[𝑖] - i, j, k = C[Li].I - V = v3D[i,j,k] + L𝑖 = Lwet[𝑖] + C𝑖 = C[L𝑖] + i, j, k = C𝑖.I + V = v3D[C𝑖] a = area2D[i,j] # From Bottom - if k < nz - k′ = k + 1 - 𝑗B = Lwet3D[i,j,k′] - if !ismissing(𝑗B) && Ω[𝑗B] # only continue if inside Ω - d = abs(zt[k] - zt[k′]) - pushTmixingvalues!(𝑖s, 𝑗s, Tvals, 𝑖, 𝑗B, κV, a, d, V) - end - end + C𝑗B = k₊₁(C𝑖, gridtype) + if !isnothing(C𝑗B) + 𝑗B = Lwet3D[C𝑗B] + if !ismissing(𝑗B) && Ω[𝑗B] # only continue if inside Ω + _, _, k′ = C𝑗B.I + d = abs(zt[k] - zt[k′]) + pushTmixingvalues!(𝑖s, 𝑗s, Tvals, 𝑖, 𝑗B, κV, a, d, V) + end + end # From Top - if k > 1 - k′ = k - 1 - 𝑗T = Lwet3D[i,j,k′] - if !ismissing(𝑗T) && Ω[𝑗T] # only continue if inside Ω - d = abs(zt[k] - zt[k′]) - pushTmixingvalues!(𝑖s, 𝑗s, Tvals, 𝑖, 𝑗T, κV, a, d, V) - end - end + C𝑗T = k₋₁(C𝑖, gridtype) + if !isnothing(C𝑗T) + 𝑗T = Lwet3D[C𝑗T] + if !ismissing(𝑗T) && Ω[𝑗T] # only continue if inside Ω + _, _, k′ = C𝑗T.I + d = abs(zt[k] - zt[k′]) + pushTmixingvalues!(𝑖s, 𝑗s, Tvals, 𝑖, 𝑗T, κV, a, d, V) + end + end end return 𝑖s, 𝑗s, Tvals end \ No newline at end of file diff --git a/src/preprocessing.jl b/src/preprocessing.jl index 23d272e..faf2d15 100644 --- a/src/preprocessing.jl +++ b/src/preprocessing.jl @@ -159,7 +159,9 @@ function makemodelgrid(; areacello, volcello, lon, lat, lev, lon_vertices, lat_v 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) - return (; area2D, v3D, DZT3d, lon_vertices, lat_vertices, lon, lat, zt, edge_length_2D, distance_to_edge_2D) + gridtype = gridtopology(lon_vertices, lat_vertices, zt) + + return (; area2D, v3D, DZT3d, lon_vertices, lat_vertices, lon, lat, zt, edge_length_2D, distance_to_edge_2D, gridtype) end function makeindices(v3D) @@ -190,34 +192,6 @@ end # │ # 1 ────┘ 2 # -# Given lon_vertices and lat_vertices, find the permutation -# that sorts the vertices in that order. -function vertexpermutation(lon_vertices, lat_vertices) - # Make sure the vertices are in the right shape (4, nx, ny) - @assert size(lon_vertices, 1) == size(lat_vertices, 1) == 4 - # Take the first grid cell - 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])) - # Find the common points - common_east = intersect(Set(points), Set(points_east)) - common_noth = intersect(Set(points), Set(points_north)) - # Find the indices of the common points - idx_east = findall(in(common_east), points) - idx_north = findall(in(common_noth), points) - idx3 = only(intersect(idx_east, idx_north)) # common to all 3 cells - idx2 = only(setdiff(idx_east, idx3)) # common to (i,j) and (i+1,j) only - idx4 = only(setdiff(idx_north, idx3)) # common to (i,j) and (i,j+1) only - idx1 = only(setdiff(1:4, idx2, idx3, idx4)) # only in (i,j) - return [idx1, idx2, idx3, idx4] -end - -# View form top for vlon and vlat vertices -# 4 ────┐ 3 -# │ -# 1 ────┘ 2 function vertexindices(dir) dir == :south ? (1, 2) : dir == :east ? (2, 3) : diff --git a/src/topology.jl b/src/topology.jl new file mode 100644 index 0000000..63f6ffb --- /dev/null +++ b/src/topology.jl @@ -0,0 +1,56 @@ +abstract type AbstractGrid end +struct BipolarGrid <: AbstractGrid + nx::Int64 + ny::Int64 + nz::Int64 +end +struct TripolarGrid <: AbstractGrid + nx::Int64 + ny::Int64 + nz::Int64 +end +struct UnknownGrid <: AbstractGrid + nx::Int64 + ny::Int64 + nz::Int64 +end + +function gridtopology(lon_vertices, lat_vertices, lev) + nx = size(lon_vertices, 2) + 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] + # If all the latitudes of the "northest" vertices are 90, then it's a "regular" grid with 2 poles + if all(NPlat .== 90) + return BipolarGrid(nx,ny,nz) + # Otherwise check if the north pole is split in two + elseif NPlon == rot180(NPlon) && NPlat == rot180(NPlat) + return TripolarGrid(nx,ny,nz) + else + return UnknownGrid(nx,ny,nz) + end +end + +# Default behavior for grids (assumed bipolar) +# wrap around in the i direction (longitude) +i₊₁(C, g::AbstractGrid) = C.I[1] < g.nx ? C .+ CartesianIndex(1, 0, 0) : CartesianIndex(1, C.I[2], C.I[3]) +i₋₁(C, g::AbstractGrid) = C.I[1] > 1 ? C .- CartesianIndex(1, 0, 0) : CartesianIndex(g.nx, C.I[2], C.I[3]) +# No connection by default in the j direction (latitude) +j₊₁(C, g::AbstractGrid) = C.I[2] < g.ny ? C .+ CartesianIndex(0, 1, 0) : nothing +j₋₁(C, ::AbstractGrid) = C.I[2] > 1 ? C .- CartesianIndex(0, 1, 0) : nothing +# No connection by default in the k direction (depth) +k₊₁(C, g::AbstractGrid) = C.I[3] < g.nz ? C .+ CartesianIndex(0, 0, 1) : nothing +k₋₁(C, ::AbstractGrid) = C.I[3] > 1 ? C .- CartesianIndex(0, 0, 1) : nothing + +# Special behavior for tripolar grids +j₊₁(C, g::TripolarGrid) = C.I[2] < g.ny ? C .+ CartesianIndex(0, 1, 0) : CartesianIndex(mod1(g.nx - C.I[1] + 1, g.nx), g.ny, C.I[3]) + +# error for unknown grids +i₊₁(C, ::UnknownGrid) = error("Unknown grid type") +i₋₁(C, ::UnknownGrid) = error("Unknown grid type") +j₊₁(C, ::UnknownGrid) = error("Unknown grid type") +j₋₁(C, ::UnknownGrid) = error("Unknown grid type") +k₊₁(C, ::UnknownGrid) = error("Unknown grid type") +k₋₁(C, ::UnknownGrid) = error("Unknown grid type") diff --git a/test/localbuild.jl b/test/localbuild.jl index ae4ec48..abb217b 100644 --- a/test/localbuild.jl +++ b/test/localbuild.jl @@ -491,23 +491,56 @@ end - # plot lon lat grid + # plot lon lat grid edges fig = Figure(size=(2000, 1000)) + loncut = lon_vertices[1,1,1] + # TODO: this is currently not working for plotting cells that wrap around the dateline + # It's OK for longitudes to be in the same window for checking connections, + # but it would be good to implement a simple function to fix that for plotting. + lon_edges = [ + lon_vertices[1,:,:] lon_vertices[4,:,end] + lon_vertices[2,end,:]' lon_vertices[3,end,end] + ] + lon_edges[1:end-1,:] .= mod.(lon_edges[1:end-1,:] .- loncut, 360) .+ loncut + lon_edges[end,:] .= mod1.(lon_edges[end,:] .- loncut, 360) .+ loncut + lat_edges = [ + lat_vertices[1,:,:] lat_vertices[4,:,end] + lat_vertices[2,end,:]' lat_vertices[3,end,end] + ] + ax = Axis(fig[1,1], xlabel = "longitude (°)", ylabel = "latitude (°)") + lonnorthsouth = reduce(vcat, [[col; NaN] for col in eachcol(lon_edges)]) + latnorthsouth = reduce(vcat, [[col; NaN] for col in eachcol(lat_edges)]) + loneastwest = reduce(vcat, [[row; NaN] for row in eachrow(lon_edges)]) + lateastwest = reduce(vcat, [[row; NaN] for row in eachrow(lat_edges)]) + lines!(ax, [lonnorthsouth; NaN; loneastwest], [latnorthsouth; NaN; lateastwest], linewidth = 0.5) + # Add boarders + lines!(ax, lon_edges[1,:], lat_edges[1,:], linewidth = 2, color=:black) + lines!(ax, lon_edges[end,:], lat_edges[end,:], linewidth = 2, color=:red) + lines!(ax, lon_edges[:,1], lat_edges[:,1], linewidth = 2, color=:black) + lines!(ax, lon_edges[:,end], lat_edges[:,end], linewidth = 2, color=:green) + Label(fig[0,1], text = "Grid cell edges check $model model", tellwidth = false, fontsize = 20) + ylims!(ax, (50, 91)) + fig + outputfile = joinpath(outputdir, "grid_cell_edges_check_local_$model.png") + @info "Saving lon/lat check as image file:\n $(outputfile)" + save(outputfile, fig) - loncut = (lon[1,1] + lon[end,1]) / 2 - lon = mod.(lon .- loncut, 360) .+ loncut + + # plot lon lat grid centers + fig = Figure(size=(2000, 1000)) + loncut = (lon[1,1] + lon[end,1]) / 2 + lon2 = mod.(lon .- loncut, 360) .+ loncut ax = Axis(fig[1,1], xlabel = "longitude (°)", ylabel = "latitude (°)") - lonnorthsouth = reduce(vcat, [[col; NaN] for col in eachcol(lon)]) + lonnorthsouth = reduce(vcat, [[col; NaN] for col in eachcol(lon2)]) latnorthsouth = reduce(vcat, [[col; NaN] for col in eachcol(lat)]) lines!(ax, lonnorthsouth, latnorthsouth, linewidth = 0.5) - loneastwest = reduce(vcat, [[row; NaN] for row in eachrow(lon)]) + loneastwest = reduce(vcat, [[row; NaN] for row in eachrow(lon2)]) lateastwest = reduce(vcat, [[row; NaN] for row in eachrow(lat)]) lines!(ax, loneastwest, lateastwest, linewidth = 0.5) - - Label(fig[0,1], text = "Longitude/latitude grid check $model model", tellwidth = false, fontsize = 20) + Label(fig[0,1], text = "Grid cell centers check $model model", tellwidth = false, fontsize = 20) fig - outputfile = joinpath(outputdir, "lonlat_check_local_$model.png") + outputfile = joinpath(outputdir, "grid_cell_centers_check_local_$model.png") @info "Saving lon/lat check as image file:\n $(outputfile)" save(outputfile, fig)