Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Adding interior of a field with indices in OutputWriter #3261

Open
wants to merge 7 commits into
base: main
Choose a base branch
from

Conversation

navidcy
Copy link
Collaborator

@navidcy navidcy commented Sep 9, 2023

Attempt to close #3260.

using Oceananigans

model = HydrostaticFreeSurfaceModel(grid = RectilinearGrid(size = (5, 5, 4), x = (-1e3, 1e3), y = (-1e3, 1e3), z = (-1e3, 0)))

ηᵢ(x, y, z) = 1 * exp(-(x^2 + y^2) / 2 * (2e2)^2)

set!(model, η = ηᵢ)

simulation = Simulation(model, Δt=100, stop_time = 1000)

simulation.output_writers[:surface] = JLD2OutputWriter(model, (η=model.free_surface.η,),
                                                       schedule = TimeInterval(200),
                                                       filename = "surface",
                                                       overwrite_existing = true)
JLD2OutputWriter scheduled on TimeInterval(3.333 minutes):
├── filepath: ./surface.jld2
├── 1 outputs: η
├── array type: Array{Float64}
├── including: [:grid, :coriolis, :buoyancy, :closure]
└── max filesize: Inf YiB

(with @glwagner)

@navidcy navidcy added bug 🐞 Even a perfect program still has bugs output 💾 labels Sep 9, 2023
@glwagner
Copy link
Member

Still some test to fix

restrict_to_interior(::Colon, loc, topo, N) = interior_indices(loc, topo, N)
restrict_to_interior(::Colon, ::Nothing, topo, N) = UnitRange(1, 1)
restrict_to_interior(index::UnitRange, ::Nothing, topo, N) = UnitRange(1, 1)
restrict_to_interior(::Colon, ::Colon, loc, topo, N) = interior_indices(loc, topo, N)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

here we can probably just return Colon?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No, basically this returns either 1:N (for center fields) or 1:N+1 (for face fields) --- the range of interior points.

This convenience feature is a little annoying. Basically the user can control index ranges through two arguments. One specifies indices. The other specifies whether to chop halos or not. If the user askes for index::Colon but no halos, we have to return something like 1:N when the field covers all indices. The complication addressed in this PR is when the available indices are restricted (ie the output field is windowed)

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, (maybe not in this PR) but I wonder if we can make with_halo=true the default. An example of where this "restricing to the interior" might have unintended (or silent) consequences is when we try to compute derivatives on post-processed fields and find out that we cannot express derivatives across the whole field but we are restricted to one point less because the field has been effectively chopped off to only the interior silently (I am guilty of something like this happening to me). This seems more like a behavior we want to explicitly decide to have rather than it being implicitly enforced as a default.

@@ -73,7 +73,7 @@ function default_dimensions(output, grid, indices, with_halos)
indices = Dict(name => validate_indices(indices, locs[name], grid) for name in keys(locs))

if !with_halos
indices = Dict(name => restrict_to_interior.(indices[name], locs[name], topo, size(grid))
indices = Dict(name => restrict_to_interior.(indices[name], indices[name], locs[name], topo, size(grid))
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@glwagner I did this here as a quick workaround.

But why the NetCDFOutputWriter calls restrict_to_interior? Isn't this the job of the construct_output method to do that?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

My memory is hazy and I don't work with it much so I could make a mistake but I think NetCDFOutputWriter needs to have indices to generate dimensions (unlike JLD2OutputWriter).

restrict_to_interior(::Colon, ::Colon, ::Nothing, topo, N) = UnitRange(1, 1)
restrict_to_interior(target_index::UnitRange, available_index, ::Nothing, topo, N) = UnitRange(1, 1)
restrict_to_interior(target_index::Colon, available_index::Colon, loc, topo, N) = interior_indices(loc, topo, N)
restrict_to_interior(target_index::Colon, available_index::Colon, ::Nothing, topo, N) = UnitRange(1, 1)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Doesn't follow YASGuide but we can make an exception...

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Damn.. Sorry... I was trying to make things cleaner by having the positional argument names

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ahaha no its ok, lets keep it.

@glwagner
Copy link
Member

glwagner commented Oct 5, 2023

Does this supercede #2771 ?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug 🐞 Even a perfect program still has bugs output 💾
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Bug when adding the interior of fields with non-trivial indices in OutputWriters
3 participants