how to output the free-surface

[francispoulin]

I have modified an example to try and output the free-surface height, but unfortunately it gives the error that "Error accessing eta at index 4", which I guess occur because eta is a 2D field.

What is the proper syntax to include eta as an output variable?

[glwagner]

@simone-silvestri how do you do it?

[francispoulin]

In case it is helpful to see what we tried, and did not work, here is our attempt that follows the examples.

outputs = merge(model.velocities, model.tracers, (η=η, ζ=ζ, δ=δ))
simulation.output_writers[:fields] = JLD2OutputWriter(model, outputs;                                                      filename = prefix * "_fields",                                                      schedule = TimeInterval(24hours),                                                      overwrite_existing = true)

[glwagner]

What are the other arguments to JLD2OutputWriter? What is the full error?

[francispoulin]

The arguments are pretty standard and do work in the absence of the η=η term in the first line.

The error is the following,

ERROR: BoundsError: attempt to access 16×36×1 Array{Float64, 3} at index [4:13, 4:33, 4:34]

[simone-silvestri]

The free surface lives only at grid.Nz+1 and the output writer is trying to take the data of the free surface in the interior (I guess your grid is 10x30x30).

In general output writers do not work at the moment with a combination of windowed fields and full fields (I think I wrote an issue somewhere about that). This is something we have to solve, it should not be like that. But until this problem is solved, you can write an output writer for the free surface alone

[simone-silvestri]

Actually, I just noticed that the issue lies in a function that restricts the output to the interior. Since that function is called only when you want to output without halos, if you pass with_halo = true to the OutputWriter it should also work with windowed fields

[glwagner]

Are you using with_halos = false @francispoulin ? (That's why I was hoping we could see the code, since both true or false are "standard")

[francispoulin]

Here is a minumum working example. Sorry it took me a while to post it.

using Oceananigans
using NCDatasets

grid = RectilinearGrid(CPU(),
    size=(4, 4, 4), halo=(3, 3, 3),
    x = (0, 1), y = (0, 1), z = (-1, 0),
    topology=(Periodic, Periodic, Bounded)
)

model = HydrostaticFreeSurfaceModel(; grid)
simulation = Simulation(model, Δt=0.1, stop_time=0.2)

outputs = merge(model.velocities, (eta=model.free_surface.η,))

simulation.output_writers[:fields] = NetCDFOutputWriter(
    model, outputs;
    filename = "tmp",
    schedule = TimeInterval(1),
    with_halos = true,
    overwrite_existing = true
)

run!(simulation)

When I try with_halos = true, I get an error that starts off with this,

ERROR: LoadError: NetCDF error: size mismatch for variable 'eta' in file './tmp.nc'. Trying to write (10, 10, 1) elements while [10, 10, 11, 1] are expected (NetCDF error code: -57)

When I try `with-halos=false', I get the error that I mentioned before,

ERROR: LoadError: BoundsError: attempt to access 10×10×1 Array{Float64, 3} at index [4:7, 4:7, 4:8]

Any advice on how to fix this?

[glwagner]

Nice thanks for the very clear MWE. I think this is a bug that the NetCDFOutputWriter incorrectly infers the indices / sizes for windowed fields like η. Probably we should make this an issue.

As a stop gap, you can use a function and manually provide the dimensions, eg something like

η_interior(model) = interior(model.free_surface.η, :, :, model.grid.Nz+1)

# This might also work:
η_interior(model) = interior(model.free_surface.η)

Then you have to pass the dimensions to NetCDFOutputWriter according to its docstring for function output...

[francispoulin]

Thanks! I will give that a try later today and let you know how it goes.

[francispoulin]

I tried a similar example using JLD2 and it worked. See below. Since I use this much more then NeteCDF, I am happy to contiue doing this. I will close this and we can revisit how eta is saved in a future PR I suppose.

using Oceananigans

grid = RectilinearGrid(CPU();   
                           size = (4, 4, 4),
                           halo = (3, 3, 3),
                              x = (0,  1),
                              y = (0,  1),
                              z = (-1, 1),
                       topology = (Periodic, Periodic, Bounded))

model = HydrostaticFreeSurfaceModel(grid = grid)
simulation = Simulation(model, Δt=0.1, stop_time=0.5)

prefix = "tmp"
outputs = merge(model.velocities, (η=model.free_surface.η,))

simulation.output_writers[:fields] = JLD2OutputWriter(model, outputs;
                                                      filename = prefix * "_fields",
                                                      schedule = TimeInterval(24hours),
                                                      overwrite_existing = true,
                                                      with_halos=true)

run!(simulation)

[NoraLoose]

Has anyone been able to output the free surface using the NetCDFOutputWriter? For context:

• I’m using ClimaOcean on top of Oceananigans.
• I don't necessarily want to output the halos, but can live with them if that makes things easier.

Following the suggestions above, I tried the following:

free_surface = (eta=ocean.model.free_surface.η,)

ocean.output_writers[:average_free_surface_writer] = NetCDFOutputWriter(ocean.model, free_surface;
                                                  schedule = AveragedTimeInterval(1days; window=1days, stride=1),
                                                  filename = "average_free_surface",
                                                  with_halos = true,
                                                  overwrite_existing = true,
                                                  array_type = Array{Float32})

The lines above execute without errors - they return:

NetCDFOutputWriter scheduled on TimeInterval(1 day):
├── filepath: average_free_surface.nc
├── dimensions: zC(54), zF(55), xC(374), yF(165), xF(374), yC(164), time(0)
├── 1 outputs: eta averaged on AveragedTimeInterval(window=1 day, stride=1, interval=1 day)
└── array type: Array{Float32}
├── file_splitting: NoFileSplitting
└── file size: 21.2 KiB

However, when I run the simulation

run!(simulation)

I encounter the following error:

DimensionMismatch: new dimensions (374, 164, 55, 1) must be consistent with array size 61336

I can provide the full stack trace if that would help.

I'm interested in trying something like this:

As a stop gap, you can use a function and manually provide the dimensions, eg something like

η_interior(model) = interior(model.free_surface.η, :, :, model.grid.Nz+1)

Is there a simple way to adapt this function when using ClimaOcean, where model is now ocean.model (where the latter is certainly not a valid function argument)?

[glwagner]

Is there a simple way to adapt this function when using ClimaOcean, where model is now ocean.model (where the latter is certainly not a valid function argument)?

If η_interior(model) doesn't work, then your model may be an OceanSeaIceModel in which case you should use

η_interior(cm) = interior(cm.ocean.model.free_surface.η, :, :, cm.ocean.model.grid.Nz+1)