FieldTimeSeries support for NetCDF

[tomchor]

@ali-ramadhan and I are interesting in expanding FieldTimeSeries to support NetCDF.

None of us is super familiar with FTS, but it seems like, for the most part, all we need would be to figure out a way to reconstruct grids based on NetCDF output and figure out how to deal with memory access if the data is PartlyInMemory. Is that correct or are there other items that we pay attention to?

For the former we can output the relevant grid parameters to file and reconstruct it in a way that's already done for JLD2 whenever the grid can't be restored directly:

Oceananigans.jl/src/OutputReaders/field_time_series.jl

Lines 489 to 529 in 82ad840

	if grid isa RectilinearGrid # we can try...
	@info "Initial attempt to transfer grid to $architecture failed."
	@info "Attempting to reconstruct RectilinearGrid on $architecture manually..."
	Nx = file["grid/Nx"]
	Ny = file["grid/Ny"]
	Nz = file["grid/Nz"]
	Hx = file["grid/Hx"]
	Hy = file["grid/Hy"]
	Hz = file["grid/Hz"]
	xᶠᵃᵃ = file["grid/xᶠᵃᵃ"]
	yᵃᶠᵃ = file["grid/yᵃᶠᵃ"]
	zᵃᵃᶠ = file["grid/zᵃᵃᶠ"]
	x = file["grid/Δxᶠᵃᵃ"] isa Number ? (xᶠᵃᵃ[1], xᶠᵃᵃ[Nx+1]) : xᶠᵃᵃ
	y = file["grid/Δyᵃᶠᵃ"] isa Number ? (yᵃᶠᵃ[1], yᵃᶠᵃ[Ny+1]) : yᵃᶠᵃ
	z = file["grid/Δzᵃᵃᶠ"] isa Number ? (zᵃᵃᶠ[1], zᵃᵃᶠ[Nz+1]) : zᵃᵃᶠ
	topo = topology(grid)
	N = (Nx, Ny, Nz)
	# Reduce for Flat dimensions
	domain = Dict()
	for (i, ξ) in enumerate((x, y, z))
	if topo[i] !== Flat
	if !(ξ isa Tuple)
	chopped_ξ = ξ[1:N[i]+1]
	else
	chopped_ξ = ξ
	end
	sξ = (:x, :y, :z)[i]
	domain[sξ] = chopped_ξ
	end
	end
	size = Tuple(N[i] for i=1:3 if topo[i] !== Flat)
	halo = Tuple((Hx, Hy, Hz)[i] for i=1:3 if topo[i] !== Flat)
	RectilinearGrid(architecture; size, halo, topology=topo, domain...)
	else
	throw(err)
	end

We were also wondering if this is worth doing over one PR and a few small ones.

cc @glwagner

[glwagner]

For a very first PR, it is only necessary to figure out how to reconstruct the grid and location. We need the grid and location in order to do any meaningful post-processing, such as computing derivatives, integrals, etc, with the saved data. If indices are not supported, we can still build 3D data. Making it useful in the long term will probably require supporting indices though. Boundary conditions are less important, but that can be put on a longer term TODO list.

I would focus completely on one backend --- probably InMemory. We may want to refactor the backends a little bit. For example, partly in memory should be support automatically if we do this right.

[glwagner]

Small PRs are always better because they lead to better design and easier review.

[ali-ramadhan]

Maybe a good starting point would be to add enough functionality so that a RectilinearGrid can be reconstructed? And tests for a the reconstruction of a simple FieldTimeSeries. This could be checked pretty thoroughly with some unit tests.

Small PRs are always better because they lead to better design and easier review.

Definitely agree! I might instead say that a PR should have one clear purpose. Perhaps the narrower the better. I suppose even narrowly-scoped PRs can grow quite large though.

[glwagner]

Maybe a good starting point would be to add enough functionality so that a RectilinearGrid can be reconstructed?

That works! But don't reinvent the wheel or implement something too specific. I would start by looking at the existing functionality for reconstructing grids eg

Oceananigans.jl/src/Grids/rectilinear_grid.jl

Lines 374 to 394 in 82ad840

	function constructor_arguments(grid::RectilinearGrid)
	arch = architecture(grid)
	FT = eltype(grid)
	args = Dict(:architecture => arch, :number_type => eltype(grid))
	# Kwargs
	topo = topology(grid)
	size = (grid.Nx, grid.Ny, grid.Nz)
	halo = (grid.Hx, grid.Hy, grid.Hz)
	size = pop_flat_elements(size, topo)
	halo = pop_flat_elements(halo, topo)
	kwargs = Dict(:size => size,
	:halo => halo,
	:x => cpu_face_constructor_x(grid),
	:y => cpu_face_constructor_y(grid),
	:z => cpu_face_constructor_z(grid),
	:topology => topo)
	return args, kwargs
	end

I think you have to save the topology and location as strings. Does NetCDF support tuples?

[glwagner]

You could have constructor_arguments(filename::String) and then something like

function reconstruct_grid(filename)
    args, kwargs = constructor_arguments(filename)
    GridType = get_grid_type(filename)
    arch = args[:architecture]
    FT = args[:number_type]
    return GridType(arch, FT; kwargs...)
end

I think that will work for lat lon and rectilinear?

For OrthogonalSphericalShellGrid you probably just have to save all metrics, not sure a way around that.

[glwagner]

This is cleaner / nicer than saving all the metrics down and loading them back, because you might end up with loaded grids that don't match the original (eg ranges get converted to arrays?)