How to output the grid metrics?

[NoraLoose]

Is there a way to write the grid metrics (distances, areas, volumes) to file?

[glwagner]

@tomchor and @ali-ramadhan may be able to answer this for NetCDFOutputWriter. I think some metrics are saved but not all because there's been talk of improving the situation.

For JLD2OutputWriter you might want to write an init function. However, if you're staying in julia it is simpler to just rebuild the grid, eg writing something like

ut = FieldTimeSeries(filename, "u") # assumes you saved a field called "u"
grid = ut.grid
Δx = xspacings(grid, Face(), Center(), Center())

I think there is an issue somewhere about compatibility with xarray / xgcm if that is specifically what you are looking for.

[glwagner]

I believe the grid is serialized in the Checkpointer but you may have to open the file with JLD2 to get it. I can figure that out.

[glwagner]

This should work:

julia> using JLD2, Oceananigans

julia> filename = "north_atlantic_calibration_1.0_i0_e0_iteration97200.jld2"
"north_atlantic_calibration_1.0_i0_e0_iteration97200.jld2"

julia> file = jldopen(filename)
JLDFile /Users/gregorywagner/Projects/ClimaOceanCalibration.jl/test/north_atlantic_calibration_1.0_i0_e0_iteration97200.jld2 (read-only)
 ├─🔢 checkpointed_properties
 ├─🔢 grid
 ├─🔢 particles
 ├─🔢 clock
 ├─🔢 coriolis
 ├─🔢 buoyancy
 ├─🔢 closure
 ├─📂 timestepper
 │  └─ ⋯ (2 more entries)
 └─ ⋯ (9 more entries)

julia> grid = file["grid"]
100×60×60 ImmersedBoundaryGrid{Float64, Bounded, Bounded, Bounded} on CPU with 7×7×3 halo:
├── immersed_boundary: GridFittedBottom(mean(z)=-2259.77, min(z)=-6000.0, max(z)=0.0)
├── underlying_grid: 100×60×60 LatitudeLongitudeGrid{Float64, Bounded, Bounded, Bounded} on CPU with 7×7×3 halo and with precomputed metrics
├── longitude: Bounded  λ ∈ [-100.0, 0.0]  regularly spaced with Δλ=1.0
├── latitude:  Bounded  φ ∈ [0.0, 60.0]    regularly spaced with Δφ=1.0
└── z:         Bounded  z ∈ [-6000.0, 0.0] variably spaced with min(Δz)=5.24961, max(Δz)=438.409

julia> close(file)

A key caveat is that the version of Oceananigans may need to be the same. Specifically grids were very recently updated between 0.94 and 0.95 to support z* vertical coordinate. So you won't be able to deserialize a grid that was saved pre-0.94 if your project is updated to version 0.95.

[glwagner]

The grids do not change very often otherwise...

[NoraLoose]

For JLD2OutputWriter you might want to write an init function. However, if you're staying in julia it is simpler to just rebuild the grid, eg writing something like

ut = FieldTimeSeries(filename, "u") # assumes you saved a field called "u"
grid = ut.grid
Δx = xspacings(grid, Face(), Center(), Center())

I guess there is no easy equivalent for getting at the grid area, is there? Specifically, I have a LatitudeLongitudeGrid, which I rebuild from the pickup file as you describe above, but then I would still like to convert

grid.Azᶠᶠᵃ

to a Field{Center, Center, Center}(grid) field.

[simone-silvestri]

We don't have those utilities yer at the moment.
You can still build the metrics using a KernelFunctionOperation as such:

ut = FieldTimeSeries(filename, "u") # assumes you saved a field called "u"
grid = ut.grid
Az = Field(KernelFunctionOperation{Face, Face, Center}(zarea, grid, Face(), Face(), Center()))
compute!(Az)

[tomchor]

Like @glwagner mentioned, we don't currently output all the grid metrics to NeTCDF, but @ali-ramadhan and I are planning on dedicating some time to implementing this soon-ish. We've been finding it hard to find the time though, se let me know if you're interested in helping out :)

In the meantime, this is what I've been personally using whenever I need xgcm integration: https://pastebin.com/vSbxz7mU
I don't remember if I adapted this to #3974 (if any adaptation is needed), so it might need some small tweaks.

[NoraLoose]

No, xgcm does not look for the name volcello. Usually, you tell xgcm what the names of your grid metrics are, see cell 4 in this notebook. There are some checks within xgcm to automatically figure out which of the variables in your dataset are the metrics, so this is how it might gotten to volcello. But there is nothing special about this name, so we don't need to adopt it (luckily)!

[glwagner]

I was trying to figure out if volcello is some kind of standard name defined by CF or S-grid but it doesn't look like it is.

[NoraLoose]

Haha, looking back at the mysterious volcello-notebook it seems to me that there is no computation in that notebook that actually needs the grid cell volume. I don't think volcello is part of the Grid object that is created in cell 4.

[tomchor]

I understand your point @tomchor. But if only half of the cell in question is wet and the other half is dry, then this point is most likely masked in the underlying fluid model, and the value at that point would not enter the average that you compute in your example.

So both solutions (leaving boundary distances as they are and dividing boundary distances by two) will lead to the same result in your computation. This is why I would go with the simpler (less confusing) way to leave the boundary distances untouched. I think this problem should rather be addressed by providing the masks.

I see your point. So the mask would be 1 in all center points, half at the edges (for face fields in a bounded topology) and zero outside the wet domain?

[glwagner]

Plus on the staggered grid I don't think it's quite right to think of the "peripheral cells" (as we have termed them in Oceananigans) as "half wet, half dry". Instead these cells either have boundary conditions (like the boundary normal velocities), or have to be diagnosed from other quantities correctly by applying the boundary conditions on those other quantities (for example the gradient of a field normal to a boundary). It's not correct to divide the metrics by two, I don't think, either (in a finite volume sense).

In other words the C-grid provides a set of rules that we can apply to understand the behavior of the peripheral points, but in no situation is it correct to think of the peripheral volumes as having half of the full cell volume.