Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add inbounds, loop along columns, fix docstring #63

Open
wants to merge 8 commits into
base: main
Choose a base branch
from
Open

Add inbounds, loop along columns, fix docstring #63

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
be4850e
add inbounds, loop along columns
navidcy Dec 20, 2024
c6c930f
fix hyperbolae/ellipses eqs to match Murray 1996
navidcy Dec 20, 2024
2ded14d
minor cleanup
navidcy Dec 20, 2024
5f29949
minor tweak
navidcy Dec 20, 2024
6978e43
some cos/sin -> cosd/sind
navidcy Dec 20, 2024
669853b
sind/cosd are custom
navidcy Dec 20, 2024
bf86f08
sind/cosd are custom
navidcy Dec 20, 2024
2e00b0d
some space
navidcy Dec 20, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
44 changes: 22 additions & 22 deletions src/generate_tripolar_coordinates.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
"""
_compute_tripolar_coordinates!(λFF, φFF, λFC, φFC, λCF, φCF, λCC, φCC,
λᶠᵃᵃ, λᶜᵃᵃ, φᵃᶠᵃ, φᵃᶜᵃ,
_compute_tripolar_coordinates!(λFF, φFF, λFC, φFC, λCF, φCF, λCC, φCC,
λᶠᵃᵃ, λᶜᵃᵃ, φᵃᶠᵃ, φᵃᶜᵃ,
first_pole_longitude,
focal_distance, Nλ)

Expand All @@ -15,39 +15,39 @@ The `focal_distance` argument is the distance from the center of the ellipses to
The family of ellipses obeys:

```math
\\frac{x²}{a² \\cosh²(ψ)} + \\frac{y²}{a² \\sinh²(ψ)} = 1
\\frac{x²}{a² \\cosh² ψ} + \\frac{y²}{a² \\sinh² ψ} = 1
```

While the family of perpendicular hyperbolae obey:
while the family of perpendicular hyperbolae obey:

```math
\\frac{x²}{a² \\cosh²(λ)} + \\frac{y²}{a² \\sinh²(λ)} = 1
\\frac{x²}{a² \\mathrm{sind}² λ} - \\frac{y²}{a² \\mathrm{cosd}² λ} = 1
```

Where ``a`` is the `focal_distance` to the center, ``λ`` is the longitudinal angle,
Above, ``a`` is the `focal_distance` to the center, ``λ`` is the longitudinal angle,
and ``ψ`` is the "isometric latitude", defined by Murray (1996) and satisfying:

```math
a \\sinh(ψ) = \\mathrm{tand}[(90 - φ) / 2]
a \\sinh ψ = \\mathrm{tand}[(90 - φ) / 2]
```

The final ``(x, y)`` points that define the stereographic projection of the tripolar
coordinates are given by:

```math
\\begin{align}
x & = a \\sinh ψ \\cos λ \\\\
y & = a \\sinh ψ \\sin λ
\\end{align}
\\begin{align*}
x & = a \\sinh ψ \\, \\mathrm{cosd} λ \\\\
y & = a \\sinh ψ \\, \\mathrm{sind} λ
\\end{align*}
```

for which it is possible to retrieve the longitude and latitude by:
from which we can recover the longitude and latitude as:

```math
\\begin{align}
λ &= - \\frac{180}{π} \\mathrm{atan}(y / x) \\\\
φ &= 90 - \\frac{360}{π} \\mathrm{atan} \\sqrt{x² + y²}
\\end{align}
\\begin{align*}
λ & = - \\frac{180}{π} \\mathrm{atan}(y / x) \\\\
φ & = 90 - \\frac{360}{π} \\mathrm{atan} \\sqrt{x² + y²}
\\end{align*}
```
"""
@kernel function _compute_tripolar_coordinates!(λFF, φFF, λFC, φFC, λCF, φCF, λCC, φCC,
Expand All @@ -72,16 +72,16 @@ for which it is possible to retrieve the longitude and latitude by:
# This makes sense, what is the longitude of the north pole? Could be anything!
# so we choose a value that is continuous with the surrounding points.
on_the_north_pole = (x == 0) & (y == 0)
north_pole_value = ifelse(i == 1, -90, 90)
north_pole_value = ifelse(i == 1, -90, 90)

λ2D[i, j] = ifelse(on_the_north_pole, north_pole_value, - 180 / π * atan(y / x))
φ2D[i, j] = 90 - 360 / π * atan(sqrt(y^2 + x^2)) # The latitude will be in the range [-90, 90]
φ2D[i, j] = 90 - 360 / π * atan(sqrt(y^2 + x^2)) # The latitude is in the range [-90, 90]

# Shift longitude to the range [-180, 180], the
# the north singularities will be located at -90 and 90
λ2D[i, j] += ifelse(i ≤ Nλ÷2, -90, 90)
# Shift longitude to the range [-180, 180],
# the north singularities are located at -90 and 90.
λ2D[i, j] += ifelse(i ≤ Nλ÷2, -90, 90)

# Make sure the singularities are at longitude we want them to be at.
# Make sure the singularities are at longitudes we want them to be at.
# (`first_pole_longitude` and `first_pole_longitude` + 180)
λ2D[i, j] += first_pole_longitude + 90
λ2D[i, j] = convert_to_0_360(λ2D[i, j])
Expand Down
88 changes: 48 additions & 40 deletions src/tripolar_grid.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,8 @@
""" a structure to represent a tripolar grid on a spherical shell """
"""
struct Tripolar{N, F, S}

A structure to represent a tripolar grid on a spherical shell.
"""
struct Tripolar{N, F, S}
north_poles_latitude :: N
first_pole_longitude :: F
Expand All @@ -22,12 +26,9 @@ const TripolarGrid{FT, TX, TY, TZ, CZ, A, Arch} = OrthogonalSphericalShellGrid{F
north_poles_latitude = 55,
first_pole_longitude = 70)

Construct a tripolar grid on a spherical shell.

!!! warning "Longitude coordinate must have even number of cells"
`size` is a 3-tuple of the grid size in longitude, latitude, and vertical directions.
Due to requirements of the folding at the north edge of the domain, the longitude size
of the grid (i.e., the first component of `size`) _must_ be an even number!
Return an `OrthogonalSphericalShellGrid` tripolar grid on the sphere. The
tripolar grid replaces the North pole singularity with two other singularities
at `north_poles_latitude` that is _less_ than 90ᵒ.

Positional Arguments
====================
Expand All @@ -47,30 +48,31 @@ Keyword Arguments
The second singularity is located at `first_pole_longitude + 180ᵒ`.
- `north_poles_latitude`: The latitude of the "north" singularities.

Return
======

An `OrthogonalSphericalShellGrid` object representing a tripolar grid on the sphere.
The north singularities are located at
!!! warning "Longitude coordinate must have even number of cells"
`size` is a 3-tuple of the grid size in longitude, latitude, and vertical directions.
Due to requirements of the folding at the north edge of the domain, the longitude size
of the grid (i.e., the first component of `size`) _must_ be an even number!

`i = 1, j = Nφ` and `i = Nλ ÷ 2 + 1, j = Nλ`
!!! info "North pole singularities"
The north singularities are located at: `i = 1`, `j = Nφ` and `i = Nλ ÷ 2 + 1`, `j = Nφ`.
"""
function TripolarGrid(arch = CPU(), FT::DataType = Float64;
size,
southernmost_latitude = -80, # The southermost `Center` latitude of the grid
southernmost_latitude = -80,
halo = (4, 4, 4),
radius = R_Earth,
z = (0, 1),
north_poles_latitude = 55,
first_pole_longitude = 70) # The second pole is at `λ = first_pole_longitude + 180ᵒ`
first_pole_longitude = 70) # second pole is at longitude `first_pole_longitude + 180ᵒ`

# TODO: change a couple of allocations here and there to be able
# TODO: Change a couple of allocations here and there to be able
# to construct the grid on the GPU. This is not a huge problem as
# grid generation is quite fast, but it might become for sub-kilometer grids
# grid generation is quite fast, but it might become slow for
# sub-kilometer resolution grids.

latitude = (southernmost_latitude, 90)
longitude = (-180, 180)

focal_distance = tand((90 - north_poles_latitude) / 2)

Nλ, Nφ, Nz = size
Expand All @@ -80,7 +82,7 @@ function TripolarGrid(arch = CPU(), FT::DataType = Float64;
throw(ArgumentError("The number of cells in the longitude dimension should be even!"))
end

# the λ and Z coordinate is the same as for the other grids,
# the λ and z coordinate is the same as for the other grids,
# but for the φ coordinate we need to remove one point at the north
# because the the north pole is a `Center`point, not on `Face` point...
topology = (Periodic, RightConnected, Bounded)
Expand Down Expand Up @@ -108,31 +110,34 @@ function TripolarGrid(arch = CPU(), FT::DataType = Float64;
φCC = zeros(Nλ, Nφ)

loop! = _compute_tripolar_coordinates!(device(CPU()), (16, 16), (Nλ, Nφ))

loop!(λFF, φFF, λFC, φFC, λCF, φCF, λCC, φCC,
λᶠᵃᵃ, λᶜᵃᵃ, φᵃᶠᵃ, φᵃᶜᵃ,
first_pole_longitude,
focal_distance, Nλ)

# We need to circshift eveything to have the first pole at the beginning of the
# We need to circshift everything to have the first pole at the beginning of the
# grid and the second pole in the middle
shift = Nλ÷4

λFF = circshift(λFF, (shift, 0))
φFF = circshift(φFF, (shift, 0))
λFC = circshift(λFC, (shift, 0))
φFC = circshift(φFC, (shift, 0))
λCF = circshift(λCF, (shift, 0))
φCF = circshift(φCF, (shift, 0))
λCC = circshift(λCC, (shift, 0))
φFF = circshift(φFF, (shift, 0))
λFC = circshift(λFC, (shift, 0))
φFC = circshift(φFC, (shift, 0))
λCF = circshift(λCF, (shift, 0))
φCF = circshift(φCF, (shift, 0))
λCC = circshift(λCC, (shift, 0))
φCC = circshift(φCC, (shift, 0))

Nx = Nλ
Ny = Nφ

# return λFF, φFF, λFC, φFC, λCF, φCF, λCC, φCC
# Helper grid to fill halo
grid = RectilinearGrid(; size = (Nx, Ny), halo = (Hλ, Hφ), topology = (Periodic, RightConnected, Flat), x = (0, 1), y = (0, 1))
grid = RectilinearGrid(; size = (Nx, Ny),
halo = (Hλ, Hφ),
x = (0, 1), y = (0, 1),
topology = (Periodic, RightConnected, Flat))
z_faces = cpu_face_constructor_z(grid)

# Boundary conditions to fill halos of the coordinate and metric terms
Expand All @@ -151,7 +156,7 @@ function TripolarGrid(arch = CPU(), FT::DataType = Float64;

lFC = Field((Face, Center, Center), grid; boundary_conditions = default_boundary_conditions)
pFC = Field((Face, Center, Center), grid; boundary_conditions = default_boundary_conditions)

lCF = Field((Center, Face, Center), grid; boundary_conditions = default_boundary_conditions)
pCF = Field((Center, Face, Center), grid; boundary_conditions = default_boundary_conditions)

Expand All @@ -174,7 +179,7 @@ function TripolarGrid(arch = CPU(), FT::DataType = Float64;
fill_halo_regions!(lCF)
fill_halo_regions!(lFC)
fill_halo_regions!(lCC)

fill_halo_regions!(pFF)
fill_halo_regions!(pCF)
fill_halo_regions!(pFC)
Expand Down Expand Up @@ -220,7 +225,7 @@ function TripolarGrid(arch = CPU(), FT::DataType = Float64;
λᶠᶜᵃ, λᶜᶜᵃ, λᶜᶠᵃ, λᶠᶠᵃ,
φᶠᶜᵃ, φᶜᶜᵃ, φᶜᶠᵃ, φᶠᶠᵃ,
radius)

# Metrics fields to fill halos
FF = Field((Face, Face, Center), grid; boundary_conditions = default_boundary_conditions)
FC = Field((Face, Center, Center), grid; boundary_conditions = default_boundary_conditions)
Expand Down Expand Up @@ -283,7 +288,7 @@ function TripolarGrid(arch = CPU(), FT::DataType = Float64;
continue_south!(Δxᶠᶜᵃ, latitude_longitude_grid.Δxᶠᶜᵃ)
continue_south!(Δxᶜᶠᵃ, latitude_longitude_grid.Δxᶜᶠᵃ)
continue_south!(Δxᶜᶜᵃ, latitude_longitude_grid.Δxᶜᶜᵃ)

continue_south!(Δyᶠᶠᵃ, latitude_longitude_grid.Δyᶠᶜᵃ)
continue_south!(Δyᶠᶜᵃ, latitude_longitude_grid.Δyᶠᶜᵃ)
continue_south!(Δyᶜᶠᵃ, latitude_longitude_grid.Δyᶜᶠᵃ)
Expand Down Expand Up @@ -323,16 +328,17 @@ function TripolarGrid(arch = CPU(), FT::DataType = Float64;
on_architecture(arch, map(FT, Azᶠᶠᵃ)),
convert(FT, radius),
Tripolar(north_poles_latitude, first_pole_longitude, southernmost_latitude))

return grid
end

# Continue the metrics to the south with LatitudeLongitudeGrid metrics
function continue_south!(new_metric, lat_lon_metric::Number)
Hx, Hy = new_metric.offsets
Nx, Ny = size(new_metric)
for i in Hx+1:Nx+Hx, j in Hy+1:1
new_metric[i, j] = lat_lon_metric

for j in Hy+1:1, i in Hx+1:Nx+Hx
@inbounds new_metric[i, j] = lat_lon_metric
end

return nothing
Expand All @@ -342,8 +348,9 @@ end
function continue_south!(new_metric, lat_lon_metric::AbstractArray{<:Any, 1})
Hx, Hy = new_metric.offsets
Nx, Ny = size(new_metric)
for i in Hx+1:Nx+Hx, j in Hy+1:1
new_metric[i, j] = lat_lon_metric[j]

for j in Hy+1:1, i in Hx+1:Nx+Hx
@inbounds new_metric[i, j] = lat_lon_metric[j]
end

return nothing
Expand All @@ -353,8 +360,9 @@ end
function continue_south!(new_metric, lat_lon_metric::AbstractArray{<:Any, 2})
Hx, Hy = - new_metric.offsets
Nx, Ny = size(new_metric)
for i in Hx+1:Nx+Hx, j in Hy+1:1
new_metric[i, j] = lat_lon_metric[i, j]

for j in Hy+1:1, i in Hx+1:Nx+Hx
@inbounds new_metric[i, j] = lat_lon_metric[i, j]
end

return nothing
Expand Down
5 changes: 3 additions & 2 deletions src/zipper_boundary_condition.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -45,8 +45,9 @@ and
c₂ == c₃
```

This is not the case for the v-velocity (or any field on the j-faces) where the last grid point is not repeated.
Because of this redundancy, we ensure consistency by substituting the redundant part of fields Centered in x, in the last row.
This is not the case for the ``v``-velocity (or any field on the `j`-faces), where
the last grid point is not repeated. Because of this redundancy, we ensure consistency
by substituting the redundant part of fields Centered in ``x``, in the last row.
"""
ZipperBoundaryCondition(sign = 1) = BoundaryCondition(Zipper(), sign)

Expand Down
Loading