Vertical index k increases as z coordinate decreases

[glwagner]

For RegularCartesianGrid, Oceananigans.jl uses a convention in which indexing for the z coordinate in arrays is treated differently from x and y: the third index (call it k) decreases with increasing z, unlike the first and second indices which increase with increasing x and y.

This means z derivative operators look different from x and y derivative operators. For example, a z differential is defined as positive for the difference between k and k+1.

@inline function δz_e2f(f, Nz, i, j, k)
    if k == Nz
        @inbounds return f[i, j, k]
    else
        @inbounds return f[i, j, k] - f[i, j, k+1]
    end
end

from here.

I recognize that some users, especially oceanographers, appreciate this feature in their output --- in particular because k=1 corresponds to the 'surface', which has a special status in (some) oceanographic applications.

However, I believe this 'reverse indexing' convention for the third direction is potentially confusing to developers and users with backgrounds in physics or engineering and may lead to bugs, slow down development, and stymie super-users who need to write code that indexes into arrays.

I propose that we consider using the convention where indices increase with increasing coordinates in all directions (for all coordinate systems), so that for a Cartesian coordinate system z and k increase upwards, like 'x', 'y', 'i' and 'j'. I believe this behavior is probably expected for the majority of developers of fluid codes in general, and also for a large fraction of users.

To satisfy users who prefer the oceanographic indexing convention in their output, we can provide users with the option (and make it default if we wish) to flip output and grids on their head before saving, so that k=1 corresponds to the surface in output.

@ali-ramadhan @christophernhill @jm-c @SandreOuza

[ali-ramadhan]

I'll let more experienced developers (@jm-c @christophernhill ) comment on this but I do agree it would be nice if the z operators looked like the x and y operators. Whether this leads to issues when e.g. bathymetry is included I'm not sure.

[johncmarshall54]

J-M mentioned to me that he had seen mistakes made in mitgcm contributions because k decreases as z decreases. One argument for going with Greg's plan. All the operators look the same in each direction is k increases upwards. What happens in ClimaAtmos?

[johncmarshall54]

Correction: meant to say k decreases as z increases in mitgcm. We should go the other way.

[glwagner]

@lcw @jkozdon can you weigh in on how ClimaAtmos.Dycore treats its vertical coordinate (if there is a "treatment")? Is a particular index somehow reserved for the land/ocean-atmosphere boundary? (I'm not sure how this works on unstructured meshes...)

Ocean models "traditionally" use a convention in which the vertical index acts like "depth" and increases downwards (and sets the ocean surface at index 1), rather than increasing upwards in the "positive direction" like other index-coordinate pairs (therefore setting the ocean surface at, say, nz+1). We think that a common convention across all coordinates may make more sense.

[jkozdon]

Our planned treatment is that the vertical index will increase with increasing distance from the earth surface. This is only really relevant for us in the 1-D column models, since for the volume solves everything is thought of as unstructured, and any such rotations would be handled through the metric terms.

[lcw]

Good question! For our high-order accurate discontinuous polynomial elements there are at least two levels of possible structure. For example, there is the order of the elements and the order of the degrees-of-freedom on each element.

The degrees-of-freedom on each element in ClimaAtmos.Dycore are in Cartesian order where increasing index corresponds with positive direction on the reference element. Currently, the slowest running index is the vertical index. I can see an argument for changing this to the fastest index as the elements are usually aligned so that the vertical direction on the reference elements corresponds to the physical vertical direction.

As far as the elements go we just recently added a semi-structured ordering where elements are grouped in vertical stacks. This was done to keep vertical columns on the same MPI rank. Right now the positive direction is "altitude" and increases with element number. That was just an arbitrary decision and can change if needed for the ocean.

[ali-ramadhan]

Thanks for the feedback! We were also thinking that if the ocean and atmosphere shared the same vertical indexing then maybe using some shared parameterizations would be easy but as pointed out, vertical indexing isn't necessarily a rigid thing in DG. Also sounds like the ocean and atmosphere don't need to agree as they're arbitrary decisions anyways.

Interesting point about setting the vertical index to be the fastest running index to help speed up vertical integrals. We only have one vertical integral in the time stepping core so not sure if it'll make a big difference on the ocean side but something to think about! Might be relevant to #21.

We decided to switch our convention and use k=1 at the bottom. I'll submit a PR soon with this change.