From 9b227a474c6ec0d3e55b54eb07e8620b0eaf2288 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 21 Nov 2024 23:47:00 -0800 Subject: [PATCH] (0.94.3) Vector rotation for immersed boundary grids (#3939) * add this correction * add a rotation operator * add test * revert * update rotation * fix tests * actually tests are ok * is this the correct direction? * better vector rotation * better rotation * fix tests * remove test for the moment * bump minor version --- Project.toml | 2 +- .../immersed_grid_metrics.jl | 9 ++- src/Operators/vector_rotation_operators.jl | 62 ++++++++++--------- 3 files changed, 41 insertions(+), 32 deletions(-) diff --git a/Project.toml b/Project.toml index 905404df86..cec2f85661 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Oceananigans" uuid = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09" authors = ["Climate Modeling Alliance and contributors"] -version = "0.94.2" +version = "0.94.3" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" diff --git a/src/ImmersedBoundaries/immersed_grid_metrics.jl b/src/ImmersedBoundaries/immersed_grid_metrics.jl index 466ad4fa2f..48cb239768 100644 --- a/src/ImmersedBoundaries/immersed_grid_metrics.jl +++ b/src/ImmersedBoundaries/immersed_grid_metrics.jl @@ -1,7 +1,7 @@ using Oceananigans.AbstractOperations: GridMetricOperation import Oceananigans.Grids: coordinates -import Oceananigans.Operators: Δzᵃᵃᶠ, Δzᵃᵃᶜ +import Oceananigans.Operators: Δzᵃᵃᶠ, Δzᵃᵃᶜ, intrinsic_vector, extrinsic_vector # Grid metrics for ImmersedBoundaryGrid # @@ -31,3 +31,10 @@ coordinates(grid::IBG) = coordinates(grid.underlying_grid) @inline Δzᵃᵃᶠ(i, j, k, ibg::IBG) = Δzᵃᵃᶠ(i, j, k, ibg.underlying_grid) @inline Δzᵃᵃᶜ(i, j, k, ibg::IBG) = Δzᵃᵃᶜ(i, j, k, ibg.underlying_grid) + +# Extend both 2D and 3D methods +@inline intrinsic_vector(i, j, k, ibg::IBG, u, v) = intrinsic_vector(i, j, k, ibg.underlying_grid, u, v) +@inline extrinsic_vector(i, j, k, ibg::IBG, u, v) = extrinsic_vector(i, j, k, ibg.underlying_grid, u, v) + +@inline intrinsic_vector(i, j, k, ibg::IBG, u, v, w) = intrinsic_vector(i, j, k, ibg.underlying_grid, u, v, w) +@inline extrinsic_vector(i, j, k, ibg::IBG, u, v, w) = extrinsic_vector(i, j, k, ibg.underlying_grid, u, v, w) \ No newline at end of file diff --git a/src/Operators/vector_rotation_operators.jl b/src/Operators/vector_rotation_operators.jl index 7ac9da502c..f42dc70466 100644 --- a/src/Operators/vector_rotation_operators.jl +++ b/src/Operators/vector_rotation_operators.jl @@ -57,30 +57,31 @@ _intrinsic_ coordinate systems are equivalent. However, for other grids (e.g., f # 2D vectors @inline function intrinsic_vector(i, j, k, grid::OrthogonalSphericalShellGrid, uₑ, vₑ) - φᶜᶠᵃ₊ = φnode(i, j+1, 1, grid, Center(), Face(), Center()) - φᶜᶠᵃ₋ = φnode(i, j, 1, grid, Center(), Face(), Center()) - Δyᶜᶜᵃ = Δyᶜᶜᶜ(i, j, 1, grid) + φᶠᶠᵃ⁺⁺ = φnode(i+1, j+1, 1, grid, Face(), Face(), Center()) + φᶠᶠᵃ⁺⁻ = φnode(i+1, j, 1, grid, Face(), Face(), Center()) + φᶠᶠᵃ⁻⁺ = φnode(i, j+1, 1, grid, Face(), Face(), Center()) + φᶠᶠᵃ⁻⁻ = φnode(i, j, 1, grid, Face(), Face(), Center()) - # θᵢ is the rotation angle between intrinsic and extrinsic reference frame - Rcosθᵢ = deg2rad(φᶜᶠᵃ₊ - φᶜᶠᵃ₋) / Δyᶜᶜᵃ - - φᶠᶜᵃ₊ = φnode(i+1, j, 1, grid, Face(), Center(), Center()) - φᶠᶜᵃ₋ = φnode(i, j, 1, grid, Face(), Center(), Center()) - Δxᶜᶜᵃ = Δxᶜᶜᶜ(i, j, 1, grid) + Δyᶠᶜᵃ⁺ = Δyᶠᶜᶜ(i+1, j, 1, grid) + Δyᶠᶜᵃ⁻ = Δyᶠᶜᶜ(i, j, 1, grid) + Δxᶜᶠᵃ⁺ = Δxᶜᶠᶜ(i, j+1, 1, grid) + Δxᶜᶠᵃ⁻ = Δxᶜᶠᶜ(i, j, 1, grid) - Rsinθᵢ = - deg2rad(φᶠᶜᵃ₊ - φᶠᶜᵃ₋) / Δxᶜᶜᵃ + # θᵢ is the rotation angle between intrinsic and extrinsic reference frame + Rcosθ = (deg2rad(φᶠᶠᵃ⁺⁺ - φᶠᶠᵃ⁺⁻) / Δyᶠᶜᵃ⁺ + deg2rad(φᶠᶠᵃ⁻⁺ - φᶠᶠᵃ⁻⁻) / Δyᶠᶜᵃ⁻) / 2 + Rsinθ = - (deg2rad(φᶠᶠᵃ⁺⁺ - φᶠᶠᵃ⁻⁺) / Δxᶜᶠᵃ⁺ + deg2rad(φᶠᶠᵃ⁺⁻ - φᶠᶠᵃ⁻⁻) / Δxᶜᶠᵃ⁻) / 2 # Normalization for the rotation angles - Rᵢ = sqrt(Rcosθᵢ^2 + Rsinθᵢ^2) + R = sqrt(Rcosθ^2 + Rsinθ^2) u = getvalue(uₑ, i, j, k, grid) v = getvalue(vₑ, i, j, k, grid) - cosθᵢ = Rcosθᵢ / Rᵢ - sinθᵢ = Rsinθᵢ / Rᵢ + cosθ = Rcosθ / R + sinθ = Rsinθ / R - uᵢ = u * cosθᵢ + v * sinθᵢ - vᵢ = - u * sinθᵢ + v * cosθᵢ + uᵢ = u * cosθ + v * sinθ + vᵢ = - u * sinθ + v * cosθ return uᵢ, vᵢ end @@ -97,30 +98,31 @@ end # 2D vectors @inline function extrinsic_vector(i, j, k, grid::OrthogonalSphericalShellGrid, uᵢ, vᵢ) - φᶜᶠᵃ₊ = φnode(i, j+1, 1, grid, Center(), Face(), Center()) - φᶜᶠᵃ₋ = φnode(i, j, 1, grid, Center(), Face(), Center()) - Δyᶜᶜᵃ = Δyᶜᶜᶜ(i, j, 1, grid) + φᶠᶠᵃ⁺⁺ = φnode(i+1, j+1, 1, grid, Face(), Face(), Center()) + φᶠᶠᵃ⁺⁻ = φnode(i+1, j, 1, grid, Face(), Face(), Center()) + φᶠᶠᵃ⁻⁺ = φnode(i, j+1, 1, grid, Face(), Face(), Center()) + φᶠᶠᵃ⁻⁻ = φnode(i, j, 1, grid, Face(), Face(), Center()) - # θₑ is the rotation angle between intrinsic and extrinsic reference frame - Rcosθₑ = deg2rad(φᶜᶠᵃ₊ - φᶜᶠᵃ₋) / Δyᶜᶜᵃ + Δyᶠᶜᵃ⁺ = Δyᶠᶜᶜ(i+1, j, 1, grid) + Δyᶠᶜᵃ⁻ = Δyᶠᶜᶜ(i, j, 1, grid) + Δxᶜᶠᵃ⁺ = Δxᶜᶠᶜ(i, j+1, 1, grid) + Δxᶜᶠᵃ⁻ = Δxᶜᶠᶜ(i, j, 1, grid) - φᶠᶜᵃ₊ = φnode(i+1, j, 1, grid, Face(), Center(), Center()) - φᶠᶜᵃ₋ = φnode(i, j, 1, grid, Face(), Center(), Center()) - Δxᶜᶜᵃ = Δxᶜᶜᶜ(i, j, 1, grid) - - Rsinθₑ = - deg2rad(φᶠᶜᵃ₊ - φᶠᶜᵃ₋) / Δxᶜᶜᵃ + # θᵢ is the rotation angle between intrinsic and extrinsic reference frame + Rcosθ = (deg2rad(φᶠᶠᵃ⁺⁺ - φᶠᶠᵃ⁺⁻) / Δyᶠᶜᵃ⁺ + deg2rad(φᶠᶠᵃ⁻⁺ - φᶠᶠᵃ⁻⁻) / Δyᶠᶜᵃ⁻) / 2 + Rsinθ = - (deg2rad(φᶠᶠᵃ⁺⁺ - φᶠᶠᵃ⁻⁺) / Δxᶜᶠᵃ⁺ + deg2rad(φᶠᶠᵃ⁺⁻ - φᶠᶠᵃ⁻⁻) / Δxᶜᶠᵃ⁻) / 2 # Normalization for the rotation angles - Rₑ = sqrt(Rcosθₑ^2 + Rsinθₑ^2) + R = sqrt(Rcosθ^2 + Rsinθ^2) u = getvalue(uᵢ, i, j, k, grid) v = getvalue(vᵢ, i, j, k, grid) - cosθₑ = Rcosθₑ / Rₑ - sinθₑ = Rsinθₑ / Rₑ + cosθ = Rcosθ / R + sinθ = Rsinθ / R - uₑ = u * cosθₑ - v * sinθₑ - vₑ = u * sinθₑ + v * cosθₑ + uₑ = u * cosθ - v * sinθ + vₑ = u * sinθ + v * cosθ return uₑ, vₑ end