fjp/spherical coriolis #4695
There are no files selected for viewing
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,71 @@ | ||
| using Oceananigans.Grids: LatitudeLongitudeGrid, OrthogonalSphericalShellGrid, peripheral_node, φnode | ||
| using Oceananigans.Operators: Δx_qᶜᶠᶜ, Δy_qᶠᶜᶜ, Δxᶠᶜᶜ, Δyᶜᶠᶜ, hack_sind, hack_cosd | ||
| using Oceananigans.Advection: EnergyConserving, EnstrophyConserving | ||
|
There was a problem hiding this comment. This file can be deleted right? There was a problem hiding this comment. Yes, and it is now deleted. Thanks. |
||
| using Oceananigans.BoundaryConditions | ||
| using Oceananigans.Fields | ||
| using Oceananigans.AbstractOperations: KernelFunctionOperation | ||
| using Oceananigans.ImmersedBoundaries | ||
|
|
||
| """ | ||
| struct NonhydrostaticSphericalCoriolis{S, FT} <: AbstractRotation | ||
|
|
||
| A parameter object for constant rotation around a vertical axis on the sphere. | ||
|
|
||
| This is only designed to work for a LatitudeLongitudeGrid. | ||
|
|
||
| There is no attempt to conserve energy or enstrophy at this early stage. | ||
| """ | ||
| struct NonhydrostaticSphericalCoriolis{S, FT} <: AbstractRotation | ||
| rotation_rate :: FT | ||
| end | ||
|
|
||
| """ | ||
| NonhydrostaticSphericalCoriolis([FT=Float64;] | ||
| rotation_rate = Ω_Earth, | ||
|
|
||
| Return a parameter object for Coriolis forces on a sphere rotating at `rotation_rate`. | ||
|
|
||
| Keyword arguments | ||
| ================= | ||
|
|
||
| - `rotation_rate`: Sphere's rotation rate; default: [`Ω_Earth`](@ref). | ||
| """ | ||
| function NonhydrostaticSphericalCoriolis(FT::DataType=Oceananigans.defaults.FloatType; | ||
| rotation_rate = Ω_Earth, | ||
|
|
||
| return NonhydrostaticSphericalCoriolis{FT}(rotation_rate) | ||
| end | ||
|
|
||
| Adapt.adapt_structure(to, coriolis::NonhydrostaticSphericalCoriolis) = | ||
| NonhydrostaticSphericalCoriolis(Adapt.adapt(to)) | ||
|
|
||
| @inline φᶠᶠᵃ(i, j, k, grid::LatitudeLongitudeGrid) = φnode(j, grid, face) | ||
|
|
||
| @inline fᶠᶠᵃ(i, j, k, grid) = 2 * coriolis.rotation_rate * hack_sind(φᶠᶠᵃ(i, j, k, grid)) | ||
| @inline f̃ᶠᶠᵃ(i, j, k, grid) = 2 * coriolis.rotation_rate * hack_cosd(φᶠᶠᵃ(i, j, k, grid)) | ||
|
|
||
| @inline f_ℑx_vᶠᶠᵃ(i, j, k, grid, coriolis, v) = fᶠᶠᵃ(i, j, k, grid, coriolis) * ℑxᶠᵃᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, v) | ||
| @inline f̃_ℑx_wᶠᶠᵃ(i, j, k, grid, coriolis, v) = f̃ᶠᶠᵃ(i, j, k, grid, coriolis) * ℑxᶠᵃᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, w) | ||
|
|
||
| @inline f_ℑy_uᶠᶠᵃ(i, j, k, grid, coriolis, u) = fᶠᶠᵃ(i, j, k, grid, coriolis) * ℑyᵃᶠᵃ(i, j, k, grid, Δy_qᶠᶜᶜ, u) | ||
|
|
||
| @inline f̃_ℑz_uᶠᶠᵃ(i, j, k, grid, coriolis, u) = f̃ᶠᶠᵃ(i, j, k, grid, coriolis) * ℑzᵃᶠᵃ(i, j, k, grid, Δz_qᶜᶜᶠ, u) | ||
|
|
||
| @inline x_f_cross_U(i, j, k, grid, U) = | ||
| @inbounds - ℑyᵃᶜᵃ(i, j, k, grid, f_ℑx_vᶠᶠᵃ, U[2]) * Δx⁻¹ᶠᶜᶜ(i, j, k, grid) | ||
| + ℑzᵃᶜᵃ(i, j, k, grid, f̃_ℑx_wᶠᶠᵃ, U[3]) * Δx⁻¹ᶠᶜᶜ(i, j, k, grid) | ||
|
|
||
| @inline y_f_cross_U(i, j, k, grid, U) = | ||
| @inbounds + ℑxᶜᵃᵃ(i, j, k, grid, f_ℑy_uᶠᶠᵃ, U[1]) * Δy⁻¹ᶜᶠᶜ(i, j, k, grid) | ||
|
|
||
| @inline z_f_cross_U(i, j, k, grid, U) = | ||
| @inbounds - ℑxᶜᵃᵃ(i, j, k, grid, f̃_ℑz_uᶠᶠᵃ, U[1]) * Δz⁻¹ᶜᶜᶠ(i, j, k, grid) | ||
|
|
||
| function Base.show(io::IO, nonhydrostatic_spherical_coriolis::NonhydrostaticSphericalCoriolis) | ||
| rotation_rate = hydrostatic_spherical_coriolis.rotation_rate | ||
| rotation_rate_Earth = rotation_rate / Ω_Earth | ||
| rotation_rate_str = @sprintf("%.2e s⁻¹ = %.2e Ω_Earth", rotation_rate, rotation_rate_Earth) | ||
|
|
||
| return print(io, "NonhydrostaticSphericalCoriolis", '\n', | ||
| "├─ rotation rate: ", rotation_rate_str, '\n') | ||
| end | ||
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -1,25 +1,27 @@ | ||
| using Oceananigans.Grids: LatitudeLongitudeGrid, OrthogonalSphericalShellGrid, peripheral_node, φnode | ||
| using Oceananigans.Operators: Δx_qᶜᶠᶜ, Δy_qᶠᶜᶜ, Δxᶠᶜᶜ, Δyᶜᶠᶜ, hack_sind, hack_cosd | ||
| using Oceananigans.Advection: EnergyConserving, EnstrophyConserving | ||
| using Oceananigans.BoundaryConditions | ||
| using Oceananigans.Fields | ||
| using Oceananigans.AbstractOperations: KernelFunctionOperation | ||
| using Oceananigans.ImmersedBoundaries | ||
|
|
||
| """ | ||
| struct SphericalCoriolis{S, FT} <: AbstractRotation | ||
|
|
||
| A parameter object for constant rotation around a vertical axis on the sphere. | ||
| """ | ||
| struct SphericalCoriolis{H, S, FT} <: AbstractRotation | ||
| rotation_rate :: FT | ||
| scheme :: S | ||
| HydrostaticCoriolis :: H | ||
|
There was a problem hiding this comment. Will this work? Whether or not we are using the hydrostatic approximation has to be encoded in a type, which is known at compile time. There was a problem hiding this comment. I can't say that it will work. This was a suggestion that I was trying out. I think that we need a flag to specify before we define Do you have a different idea on how to do this? There was a problem hiding this comment. i propose this API: hydrostatic_coriolis = HydrostaticSphericalCoriolis()
full_coriolis = SphericalCoriolis()for implementation, we need a type such that const HydrostaticSphericalCoriolis{S, FT} = SphericalCoriolis{S, FT, <:HydrostaticFormulation} where {S, FT}note this preserves the type structure on one can introduce the concept of the "formulation" like this struct HydrostaticFormulation end
struct NonhydrostaticFormulation end
struct SphericalCoriolis{S, FT, F}
rotation_rate :: FT
scheme :: S
formulation :: F
endThere was a problem hiding this comment. Thanks @glwagner . I'll give this a try and get back to you when I have something together. |
||
| end | ||
|
|
||
| """ | ||
| SphericalCoriolis([FT=Float64;] | ||
| rotation_rate = Ω_Earth, | ||
| scheme = EnstrophyConserving()) | ||
| HydrostaticCoriolis = | ||
|
|
||
| Return a parameter object for Coriolis forces on a sphere rotating at `rotation_rate`. | ||
|
|
||
|
|
@@ -28,71 +30,60 @@ Keyword arguments | |
|
|
||
| - `rotation_rate`: Sphere's rotation rate; default: [`Ω_Earth`](@ref). | ||
| - `scheme`: Either `EnergyConserving()`, `EnstrophyConserving()`, or `EnstrophyConserving()` (default). | ||
| - 'HydrostaticCoriolis`: `HydrostaticCoriolis` or `NonhydrostaticCoriolis`. | ||
| """ | ||
|
|
||
| function SphericalCoriolis(FT::DataType=Oceananigans.defaults.FloatType; | ||
| rotation_rate = Ω_Earth, | ||
| scheme :: S = EnstrophyConserving(FT)) where S | ||
|
|
||
| return SphericalCoriolis{S, FT}(rotation_rate, scheme) | ||
| end | ||
|
|
||
| Adapt.adapt_structure(to, coriolis::SphericalCoriolis) = | ||
| SphericalCoriolis(Adapt.adapt(to, coriolis.rotation_rate), | ||
| Adapt.adapt(to, coriolis.scheme)) | ||
|
|
||
| @inline φᶠᶠᵃ(i, j, k, grid::LatitudeLongitudeGrid) = φnode(j, grid, face) | ||
| @inline φᶠᶠᵃ(i, j, k, grid::OrthogonalSphericalShellGrid) = φnode(i, j, grid, face, face) | ||
| @inline φᶠᶠᵃ(i, j, k, grid::ImmersedBoundaryGrid) = φᶠᶠᵃ(i, j, k, grid.underlying_grid) | ||
|
|
||
| @inline fᶠᶠᵃ(i, j, k, grid, coriolis::SphericalCoriolis) = | ||
| 2 * coriolis.rotation_rate * hack_sind(φᶠᶠᵃ(i, j, k, grid)) | ||
| @inline f̃ᶠᶠᵃ(i, j, k, grid, coriolis::SphericalCoriolis) = | ||
| 2 * coriolis.rotation_rate * hack_cosd(φᶠᶠᵃ(i, j, k, grid)) | ||
|
|
||
| @inline f_ℑy_uᶠᶠᵃ(i, j, k, grid, coriolis, u) = fᶠᶠᵃ(i, j, k, grid, coriolis) * ℑyᵃᶠᵃ(i, j, k, grid, Δy_qᶠᶜᶜ, u) | ||
| @inline f̃_ℑz_uᶠᶠᵃ(i, j, k, grid, coriolis::NonhydrostaticSphericalCoriolis, u) = f̃ᶠᶠᵃ(i, j, k, grid, coriolis) * ℑzᵃᶠᵃ(i, j, k, grid, Δz_qᶜᶜᶠ, u) | ||
| @inline f̃_ℑz_uᶠᶠᵃ(i, j, k, grid, coriolis::HydrostaticSphericalCoriolis, u) = zero(grid) | ||
|
|
||
| @inline f_ℑx_vᶠᶠᵃ(i, j, k, grid, coriolis, v) = fᶠᶠᵃ(i, j, k, grid, coriolis) * ℑxᶠᵃᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, v) | ||
|
|
||
| @inline f̃_ℑx_wᶠᶠᵃ(i, j, k, grid, coriolis::NonhydrostaticSphericalCoriolis, w) = f̃ᶠᶠᵃ(i, j, k, grid, coriolis) * ℑxᶠᵃᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, w) | ||
| @inline f̃_ℑx_wᶠᶠᵃ(i, j, k, grid, coriolis::HydrostaticSphericalCoriolis, w) = zero(grid) | ||
|
|
||
| @inline x_f_cross_U(i, j, k, grid, U) = | ||
| @inbounds - ℑyᵃᶜᵃ(i, j, k, grid, f_ℑx_vᶠᶠᵃ, U[2]) * Δx⁻¹ᶠᶜᶜ(i, j, k, grid) | ||
| + ℑzᵃᶜᵃ(i, j, k, grid, f̃_ℑx_wᶠᶠᵃ, U[3]) * Δx⁻¹ᶠᶜᶜ(i, j, k, grid) | ||
|
|
||
| @inline y_f_cross_U(i, j, k, grid, U) = | ||
| @inbounds + ℑxᶜᵃᵃ(i, j, k, grid, f_ℑy_uᶠᶠᵃ, U[1]) * Δy⁻¹ᶜᶠᶜ(i, j, k, grid) | ||
|
|
||
|
|
||
| @inline z_f_cross_U(i, j, k, grid, U) = | ||
| @inbounds - ℑxᶜᵃᵃ(i, j, k, grid, f̃_ℑz_uᶠᶠᵃ, U[1]) * Δz⁻¹ᶜᶜᶠ(i, j, k, grid) | ||
|
|
||
| ##### | ||
| ##### Show | ||
| ##### | ||
|
|
||
| function Base.show(io::IO, spherical_coriolis::SphericalCoriolis) | ||
| coriolis_scheme = spherical_coriolis.scheme | ||
| rotation_rate = spherical_coriolis.rotation_rate | ||
| rotation_rate_Earth = rotation_rate / Ω_Earth | ||
| rotation_rate_str = @sprintf("%.2e s⁻¹ = %.2e Ω_Earth", rotation_rate, rotation_rate_Earth) | ||
|
|
||
| return print(io, "SphericalCoriolis", '\n', | ||
| "├─ rotation rate: ", rotation_rate_str, '\n', | ||
| "└─ scheme: ", summary(coriolis_scheme)) | ||
| end | ||
Oops, something went wrong.
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
these need to be type aliases, not abstract types
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'm happy to change this but maybe best after we decide on how to swtich from one to the other.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The user will decide which one to use in their script.