-
Notifications
You must be signed in to change notification settings - Fork 191
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
Output constructor does not work with sliced fields #2770
Comments
If I try to just output the whole of the fields in the lat lon example then it fails to be able to output the free surface displacement with the same issue. |
That was the main issue we discussed in #2740 with locating the free surface at (from #2740)
|
This is a bug though. With |
The issue here seems to be with for example for a 10x10x10 grid:
So if I try to do
because
In my opinion |
Also just to clarify, while output writers accept an argument |
Only `BoundaryCondition`s as default_boundary_conditions
I'm not sure I grasp the issue. @jagoosw it might help if you include a simple example (eg some code, and the error that gets produced) that illustrates the issue.
I agree with this. Is the problem associated with trying to "remove halos" from output? Does it go away when we write In addition to problems with Here's what we have:
A few solutions:
julia> a = rand(2, 3)
2×3 Matrix{Float64}:
0.493349 0.161279 0.244427
0.288415 0.0111998 0.543645
julia> b = view(a, :, 4:4)
ERROR: BoundsError: attempt to access 2×3 Matrix{Float64} at index [1:2, 4:4]
Stacktrace:
[1] throw_boundserror(A::Matrix{Float64}, I::Tuple{Base.Slice{Base.OneTo{Int64}}, UnitRange{Int64}})
@ Base ./abstractarray.jl:651
[2] checkbounds
@ ./abstractarray.jl:616 [inlined]
[3] view(::Matrix{Float64}, ::Function, ::UnitRange{Int64})
@ Base ./subarray.jl:177
[4] top-level scope
@ REPL[36]:1
us, vs, ws, Ts = @indices (:, :, Nz) u v w T This would be equivalent to writing something like us = Field(u, indices=(:, :, Nz))
vs = Field(v, indices=(:, :, Nz))
ws = Field(w, indices=(:, :, Nz))
Ts = Field(T, indices=(:, :, Nz)) Then, if users want to output surface fields and the free surface, they should use this method rather than the |
Here's an example:
And the error I get:
|
The problem is these lines Oceananigans.jl/src/Fields/field.jl Lines 243 to 244 in b288527
(if we assume This is not the whole story though, since the output writer is trying to view the underlying data at [4:23, 4:23, 4:23]. So once we have fixed this we have to deal with the output writers which rely on Oceananigans.jl/src/OutputWriters/output_construction.jl Lines 50 to 60 in b288527
Here we have to interpolate the |
Ok! To summarize, we currently have Oceananigans.jl/src/Grids/grid_utils.jl Line 213 in b288527
which you propose to change to parent_index_range(index::UnitRange, loc, topo, halo) = 1:(last(index) - first(index) + 1) I definitely agree that the parent index range of a field with indices Note that with this change, none of the methods for
I think @navidcy discovered a problem with |
Is this still an issue? |
I've just come across this again so I think so yes. Might try and fix soon |
Hi all,
I've come across a problem a few times with the output writer where it can't set up running with
indices
specified. This presents as a bounds error coming from theoffset_data
function:Oceananigans.jl/src/Grids/new_data.jl
Lines 33 to 48 in c73ef03
An example of this occurs when you try to run the "near_global_lat_lon" example when it tries to output surface values:
Oceananigans.jl/validation/near_global_lat_lon/near_global_quarter_degree.jl
Lines 296 to 301 in c73ef03
Essentially what I think is happening is that when we output data and specify some slice that isn't the dimensions of the grid, in this example just the surface level
grid.Nz
, the data is being correctly (? not sure we actually want the halo here) selected as e.g.-Hx:Nx+Hx
in dimensions specified as:
but just whatever indices specified by the user in the other dimensions. In theoffset_data
function I thinkii
gets the halo added in all dimensions so there's a miss match.I'm not sure how to fix this cleanly?
The text was updated successfully, but these errors were encountered: