Ax, Ay and Az are negative in a Tripolar grid

[taimoorsohail]

I am trying to compute the velocity multiplied by the grid cell area (Ax, Ay or Az in Oceananigans). However, I seem to have stumbled upon a confusing feature of these Operators in that they are negative in some places in a Tripolar grid.

As a simple example:

using ClimaOcean
using Oceananigans

arch = CPU()

# ### Grid and Bathymetry
@info "Defining grid"

Nx = Integer(360*4)
Ny = Integer(180*4)
Nz = Integer(100)

@info "Defining vertical z faces"

r_faces = exponential_z_faces(; Nz, depth=5000, h=34)
z_faces = Oceananigans.MutableVerticalDiscretization(r_faces)

@info "Defining tripolar grid"

underlying_grid = TripolarGrid(arch;
                               size = (Nx, Ny, Nz),
                               z = z_faces,
                               halo = (5, 5, 4),
                               first_pole_longitude = 70,
                               north_poles_latitude = 55)

@info "Defining bottom bathymetry"

@time bottom_height = regrid_bathymetry(underlying_grid;
                                  minimum_depth = 10,
                                  interpolation_passes = 75, # 75 interpolation passes smooth the bathymetry near Florida so that the Gulf Stream is able to flow
				                  major_basins = 2)

# For this bathymetry at this horizontal resolution we need to manually open the Gibraltar strait.
view(bottom_height, 102:103, 124, 1) .= -400

@info "Defining grid"

@time grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bottom_height); active_cells_map=true)

c = XFaceField(grid)

using Oceananigans.Operators: Ax, Δz

c .= 1

compute!(Field(c * Ax))

gives

1440×720×100 Field{Face, Center, Center} on ImmersedBoundaryGrid on CPU
├── grid: 1440×720×100 ImmersedBoundaryGrid{Float64, Periodic, Oceananigans.Grids.RightConnected, Bounded} on CPU with 5×5×4 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: ZeroFlux, north: Zipper, bottom: ZeroFlux, top: ZeroFlux, immersed: Default
├── operand: BinaryOperation at (Face, Center, Center)
├── status: time=0.0
└── data: 1450×730×108 OffsetArray(::Array{Float64, 3}, -4:1445, -4:725, -3:104) with eltype Float64 with indices -4:1445×-4:725×-3:104
    └── max=4.0224e6, min=-3.57265e6, mean=1.07532e6

where the minimum is a negative value? I get the same issue with Ay and Az. Any guidance on why this is would be much appreciated.

Also, I checked and it doesn't seem to be to do with halos:

julia> maximum(interior(compute!(Field(c * Ax))))
4.022397708204129e6

julia> minimum(interior(compute!(Field(c * Ax))))
-4.022397708203645e6

[taimoorsohail]

cc @navidcy

[navidcy]

out of curiosity, have you tried other grids?

[taimoorsohail]

Yes I tried LatitudeLongitudeGrid from the near_global_ocean example and it was all positive

[navidcy]

I trimmed down the MWE a bit so we don't need ClimaOcean.

using Oceananigans
using Oceananigans.Operators: Ax, Ay, Az

Nx, Ny, Nz = 40, 30, 10

grid = TripolarGrid(size = (Nx, Ny, Nz),
                    z = (-5000, 0),
                    first_pole_longitude = 70,
                    north_poles_latitude = 55)


f = XFaceField(grid)
set!(f, 1)

f_dAx = Field(f * Ax)

compute!(f_dAx)

fig = Figure()
ax = Axis(fig[1, 1])
hm = heatmap!(ax, interior(f_dAx, :, :, 1), colorrange = extrema(f_dA), colormap = :balance)
Colorbar(fig[1, 2], hm)
fig

[navidcy]

Also

f_dAy = Field(f * Ay)

compute!(f_dAy)

fig = Figure()
ax = Axis(fig[1, 1])
hm = heatmap!(ax, interior(f_dAy, :, :, 3), colorrange = extrema(f_dA), colormap = :balance)
Colorbar(fig[1, 2], hm)
fig

and

f_dAz = Field(f * Az)

compute!(f_dAz)

fig = Figure()
ax = Axis(fig[1, 1])
hm = heatmap!(ax, interior(f_dAz, :, :, 3), colorrange = extrema(f_dA), colormap = :balance)
Colorbar(fig[1, 2], hm)
fig

[navidcy]

Btw, if

f = CenterField(grid)

everything seems reasonable.

[simone-silvestri]

Hmmm, that is because of the northern boundary condition which, for x-face fields, is assumed to need a sign switch across the fold (velocities need a sign switch across the fold).
We can check that the underlying metrics are correct (grid.Azᶠᶜᵃ for example)

To fix this problem we need to have a definition of Scalar vs Vector, where only the latter requires a sign switch.
We had a discussion some time ago of introducing a VectorField in oceananigans, maybe this is the excuse to develop it.

[navidcy]

@simone-silvestri is it because we didn't fill the halos? E.g., if I fill the halos of f after I set it I get:

using Oceananigans
using Oceananigans.Operators: Ax, Ay, Az

Nx, Ny, Nz = 40, 30, 10

grid = TripolarGrid(size = (Nx, Ny, Nz),
                    z = (-5000, 0),
                    first_pole_longitude = 70,
                    north_poles_latitude = 55)

f = XFaceField(grid)
set!(f, 1)

using Oceananigans.BoundaryConditions: fill_halo_regions!
fill_halo_regions!(f)

f_dAx = Field(f * Ax)

compute!(f_dAx)

fig = Figure()
ax = Axis(fig[1, 1])
hm = heatmap!(ax, interior(f_dAx, :, :, 3), colorrange = extrema(f_dA), colormap = :balance)
Colorbar(fig[1, 2], hm)
fig

save("tripolar_f_Ax.png", fig)


f_dAy = Field(f * Ay)

compute!(f_dAy)

fig = Figure()
ax = Axis(fig[1, 1])
hm = heatmap!(ax, interior(f_dAy, :, :, 3), colorrange = extrema(f_dA), colormap = :balance)
Colorbar(fig[1, 2], hm)
fig

save("tripolar_f_Ay.png", fig)

f_dAz = Field(f * Az)

compute!(f_dAz)

fig = Figure()
ax = Axis(fig[1, 1])
hm = heatmap!(ax, interior(f_dAz, :, :, 3), colorrange = extrema(f_dA), colormap = :balance)
Colorbar(fig[1, 2], hm)
fig

save("tripolar_f_Az.png", fig)

[simone-silvestri]

Hah, ok this is probably unexpected. I thought that the BCs are the reason for the negative values.

[simone-silvestri]

Ok if I do this

Nx, Ny, Nz = 40, 30, 10
grid = TripolarGrid(size = (Nx, Ny, Nz),
                           z = (-5000, 0),
                           first_pole_longitude = 70,
                           north_poles_latitude = 55)

f = XFaceField(grid)
set!(f, 1)
fill_halo_regions!(f)
heatmap(interior(f, :, :, 1))

as expected I get

[simone-silvestri]

And for (for example) the Ay metric, I get

f_dAy = Field(f * Ay)
compute!(f_dAy)
heatmap(interior(f_dAy, :, :, 3))

[simone-silvestri]

@navidcy I cannot reproduce your results, the negative line is always there (provided I fix the heatmap extrema because in your snipper you are using a non defined global variable)

[navidcy]

Oh darn!
I still don't get negatives with CenterField. Look:

using Oceananigans
using Oceananigans.Operators: Ax, Ay, Az
using GLMakie

Nx, Ny, Nz = 40, 30, 10

grid = TripolarGrid(size = (Nx, Ny, Nz),
                    z = (-5000, 0),
                    first_pole_longitude = 70,
                    north_poles_latitude = 55)

f = CenterField(grid)
set!(f, 1)

using Oceananigans.BoundaryConditions: fill_halo_regions!
fill_halo_regions!(f)

f_dAx = Field(f * Ax)
compute!(f_dAx)
max_abs = maximum(abs, f_dAx)
extrema_f_dA = extrema(f_dAx)

fig = Figure()
ax = Axis(fig[1, 1])
hm = heatmap!(ax, interior(f_dAx, :, :, 1), colorrange = (-max_abs, max_abs), colormap = :balance)
Label(fig[0, :], "minimum = $(extrema_f_dA[1]), maximum = $(extrema_f_dA[2])", tellwidth=false)
Colorbar(fig[1, 2], hm)
fig

save("tripolar_f_Ax.png", fig)

f_dAy = Field(f * Ay)
compute!(f_dAy)
max_abs = maximum(abs, f_dAy)
extrema_f_dA = extrema(f_dAy)

fig = Figure()
ax = Axis(fig[1, 1])
hm = heatmap!(ax, interior(f_dAy, :, :, 1), colorrange = (-max_abs, max_abs), colormap = :balance)
Label(fig[0, :], "minimum = $(extrema_f_dA[1]), maximum = $(extrema_f_dA[2])", tellwidth=false)
Colorbar(fig[1, 2], hm)
fig

save("tripolar_f_Ay.png", fig)

f_dAz = Field(f * Az)
compute!(f_dAz)
max_abs = maximum(abs, f_dAz)
extrema_f_dA = extrema(f_dAz)

fig = Figure()
ax = Axis(fig[1, 1])
hm = heatmap!(ax, interior(f_dAz, :, :, 1), colorrange = (-max_abs, max_abs), colormap = :balance)
Label(fig[0, :], "minimum = $(extrema_f_dA[1]), maximum = $(extrema_f_dA[2])", tellwidth=false)
Colorbar(fig[1, 2], hm)
fig

save("tripolar_f_Az.png", fig)

[navidcy]

With f = XFaceField(grid) I get:

as you do.

[simone-silvestri]

yeah the difference is that CenterFields are assumed to be scalar that do not require switching values, while XFaceFields are assumed to be vector components. So this behaviour is consistent with what I was expecting. You can change the boundary conditions manually if you don't want the switching of sign

I think though, the best thing to do is to write an abstraction for vector field

[glwagner]

It's better to not assume that XFaceField(grid) (with no additional annotation) is a vector component. We can add special annotation to make vector components / add assumptions specific to velocity components within model constructors.

This is also how boundary conditions work.

[glwagner]

And just to clarify, we already have a kind of vector abstraction in that we implement for other topologies than RightFolded. The main criticism could be that its conflated with the notion of "prognostic", ie we are hard-coding the C-grid when define prognostic. For example:

Oceananigans.jl/src/BoundaryConditions/field_boundary_conditions.jl

Line 21 in b37dd18

default_prognostic_bc(::Bounded, ::Face, default) = ImpenetrableBoundaryCondition()

but

Oceananigans.jl/src/BoundaryConditions/field_boundary_conditions.jl

Line 33 in b37dd18

default_auxiliary_bc(::Bounded, ::Face) = nothing

the same logic should be applied to the implementation of boundary conditions for RightFolded.

[simone-silvestri]

Yeah, that is a bit lacking as an abstraction. We leave out fields that are auxiliary but vector components (tendencies for example). Also it does not distinguish right and left boundary condition. I think This example and the needs of a right folded boundary shows the pitfalls of this design and calls for a redesign of this abstraction.

[simone-silvestri]

I was already suspecting that this piece of code was not well designed for our goals when introducing the distributed communication that does not fit within this abstraction (i.e. we need to inject boundary conditions manually).
I think the clue that we need to change this pattern is the fact that we need to do it also for a tripolar grid at the moment

[glwagner]

I think we just have to define the default so that it is not assuming a vector, and then to specify vector where needed in the model constructor. For other cases we use the word "auxiliary" for the default and "prognostic" for the vector case (because the only prognostic fields at Face are vectors). But we can just use a different word than prognostic if that's what you're saying. I don't think we need three cases?

[glwagner]

Why does fill_halo_regions! change values inside the domain? Shouldn't it only change halos?

[navidcy]

Why does fill_halo_regions! change values inside the domain? Shouldn't it only change halos?

That's a good point! I'm also confused.

[simone-silvestri]

Why does fill_halo_regions! change values inside the domain? Shouldn't it only change halos?

by construction the half line on the right hand side of the grid, for fields that are centered in y, is redundant. So it cannot really be considered as "interior" probably. At the beginning the halos were not changing this line, that would drift in time compared to the left hand side because of finite precision errors propagating. We have decided subsequently that it is best to have only one source of truth and, therefore, to substitute those values with the ones on the left hand side here CliMA/OrthogonalSphericalShellGrids.jl#58. I think this is still the best strategy numerically. The inconsistency is that it is quite difficult to enforce on distributed tripolar grids.

[glwagner]

why does the fold overlap? can we fold the grid so that the two boundaries touch rather than becoming enmeshed?