Pickup a simulation at a certain timestep and run it in a different horizontal resolution

[iuryt]

Hi all I am running a submesoscale simulation of 1km resolution and I would like to pickup at a certain timestep (e.g. 30 days of simulation), average for each 10 km and run it in a coarsen resolution.

My confort-zone answer:

1. Read the netCDF output on Python
2. Average each 10 km using xarray and save to another nc file
3. Read it on Julia and setup this as the initial conditions for the other simulation

I am looking to ways to do that using Julia (if it is easier). My idea is to have a .jl file for each simulation.
In the file for coarsen resolution simulation I would read the output from 1km simulation, average the fields and run in the new grid.

Any idea on which methods or previous examples I could explore?

[glwagner]

We don't have horizontal regridding right now but it's a highly desirable feature to have. Is generic regridding built into xarray?

The best we can do now using existing Oceananigans tools is vertical regridding, which is provided by regrid! and implemented here:

https://github.com/CliMA/Oceananigans.jl/blob/main/src/Fields/regridding_fields.jl

this feature is experimental and so undocumented.

Simple horizontal regridding could be achieved by extending that utility to work in x and y and then supporting 3D regridding by regridding each dimension in serial.

Eventually we will likely need more general regridding between grids that are only partially overlapping and possibly not coincident (eg CubedSphere to lat-lon, etc).

I'd be curious to know what tools are out there, especially in Julia. This utility:

https://github.com/ClimateGlobalChange/tempestremap

is used for climate applications and we probably will need either to use it or write something with similar functionality, and I think

https://github.com/CliMA/ClimaCoupler.jl

either uses that now or will in the future.

[iuryt]

I will take a look at this during this week, as I think that it would be a nice PR for the source code.

[iuryt]

This is not my case, but what if the target_grid is at a higher resolution and there are no points from source_field between target_grid faces?

Oceananigans.jl/src/Fields/regridding_fields.jl

Lines 101 to 116 in 00c98a7

	# Add contribution from all full cells in the integration range
	@unroll for k_src = k₋_src:k₊_src
	@inbounds target_field[i, j, k] += source_field[i_src, j_src, k_src] * Δzᶜᶜᶜ(i_src, j_src, k_src, source_grid)
	end
	zk₋_src = znode(Center(), Center(), Face(), i_src, j_src, k₋_src, source_grid)
	zk₊_src = znode(Center(), Center(), Face(), i_src, j_src, k₊_src+1, source_grid)
	# Add contribution to integral from fractional bottom part,
	# if that region is a part of the grid.
	@inbounds target_field[i, j, k] += source_field[i_src, j_src, k₋_src - 1] * (zk₋_src - z₋)
	# Add contribution to integral from fractional top part
	@inbounds target_field[i, j, k] += source_field[i_src, j_src, k₊_src] * (z₊ - zk₊_src)
	@inbounds target_field[i, j, k] /= Δzᶜᶜᶜ(i, j, k, target_grid)

Also, this code seems to be launching (looping) for each xy point.. to run this for the cell volume, I would need to change to :xyz and adapt the code to every cell, right?

Oceananigans.jl/src/Fields/regridding_fields.jl

Lines 65 to 66 in 00c98a7

	event = launch!(arch, target_grid, :xy, _regrid!, a, b, target_grid, source_grid, source_z_faces)
	wait(device(arch), event)

[glwagner]

The algorithm is supposed to work in that case (there are tests for that) --- then k₋_src = k₊_src and the value of the "fine grid" cell is simply set equal to the value of the coarse grid.

[glwagner]

Looking at the code doesn't look hard because it basically integrates field data within target-grid cells and divide by dz, which is exactly what I need but for x and y.

Right, the simplest implementation is literally just changing a few characters: z to x or y and k to i or j.

Note though to be more generally valid we may need to work a little harder. For cuvilinear validity we should use fractional volume rather than fractional spacing. Also, we need to ingest the location of the field and use the "location-agnostic" operators like Δz(i, j, k, grid, LX, LY, LZ) (or maybe this would be volume(i, j, k, grid, LX, LY, LZ) if we are using fractional volumes). In the current code we make a strong assumption that fields are at cell centers (which is irrelevant for all grids we have, since the vertical metric doesn't depend on horizontal location, but does matter on a curvilinear or stretched grid).

[glwagner]

It's still all pretty straightforward...

[glwagner]

Might be worth mentioning two other codes I recently learned about for coarse-graining and filtering geophysical data:

• FlowSieve
• GCM-filters

These tools seem oriented towards post-processing / analysis, but I don't see why they couldn't be useful, at least in principle, for forward simulations. Maybe @NoraLoose has more to add?

[christophernhill]

@iuryt you could also look at what @gaelforget has ( https://github.com/gaelforget ). Gael has some interpolation code here https://juliaclimate.github.io/MeshArrays.jl/dev/tutorials/#Geography . He is interested in collaborations to help grow collections of tools/libraries .

[glwagner]

I guess if we moved the finite volume engine of Oceananigans into its own package, it could become a common basis / community tool for any geophysics FV calculations. I think we haven’t done that because a) it’s work and B) we’d lose some the benefits of a monorepo during the “chaos phase” of development which we might still be in.