Broadcasting slices of fields doesn't work as expected #3361

navidcy · 2023-10-24T23:21:50Z

Seems that there is something strange with how broadcasting and fields work. Here, we create two fields and try to broadcast a slice of one onto the other.

julia> using Oceananigans

julia> grid = RectilinearGrid(size=(2, 3, 4), extent = (1, 1, 1))
2×3×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── Periodic x ∈ [0.0, 1.0)  regularly spaced with Δx=0.5
├── Periodic y ∈ [0.0, 1.0)  regularly spaced with Δy=0.333333
└── Bounded  z ∈ [-1.0, 0.0] regularly spaced with Δz=0.25

julia> c1 = CenterField(grid)
2×3×4 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 2×3×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
└── data: 8×9×10 OffsetArray(::Array{Float64, 3}, -2:5, -2:6, -2:7) with eltype Float64 with indices -2:5×-2:6×-2:7
    └── max=0.0, min=0.0, mean=0.0

julia> c2 = CenterField(grid)
2×3×4 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 2×3×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
└── data: 8×9×10 OffsetArray(::Array{Float64, 3}, -2:5, -2:6, -2:7) with eltype Float64 with indices -2:5×-2:6×-2:7
    └── max=0.0, min=0.0, mean=0.0

julia> c1[1, :, 1] .= c2[2, :, 3]
ERROR: DimensionMismatch: array could not be broadcast to match destination
Stacktrace:
 [1] check_broadcast_shape
   @ ./broadcast.jl:553 [inlined]
 [2] check_broadcast_axes
   @ ./broadcast.jl:556 [inlined]
 [3] instantiate
   @ ./broadcast.jl:297 [inlined]
 [4] materialize!
   @ ./broadcast.jl:884 [inlined]
 [5] materialize!
   @ ./broadcast.jl:881 [inlined]
 [6] materialize!(dest::Field{Center, Center, Center, Nothing, RectilinearGrid{Float64, Periodic, Periodic, Bounded, Float64, Float64, Float64, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, CPU}, Tuple{UnitRange{Int64}, Colon, UnitRange{Int64}}, OffsetArrays.OffsetArray{Float64, 3, SubArray{Float64, 3, Array{Float64, 3}, Tuple{UnitRange{Int64}, Base.Slice{Base.OneTo{Int64}}, UnitRange{Int64}}, false}}, Float64, FieldBoundaryConditions{Nothing, Nothing, BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Nothing, Nothing, BoundaryCondition{Oceananigans.BoundaryConditions.Flux, Nothing}}, Nothing, Oceananigans.Fields.FieldBoundaryBuffers{Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}}, bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Nothing, typeof(identity), Tuple{OffsetArrays.OffsetVector{Float64, Vector{Float64}}}})
   @ Oceananigans.Fields ~/Research/OC7.jl/src/Fields/broadcasting_abstract_fields.jl:36
 [7] top-level scope
   @ REPL[5]:1

While the size of the slices are exactly the same, we we get an error that ERROR: DimensionMismatch: array could not be broadcast to match destination.

If instead we try to broadcast to the parent arrays then things work as expected.

julia> parent(c1[1, :, 1]) .= parent(c2[2, :, 3])
9-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0

(found with @siddharthabishnu)

The text was updated successfully, but these errors were encountered:

simone-silvestri · 2023-10-25T01:30:35Z

It has to do with view(field, i, j, k) that returns a WindowedField with 3 dimensions while getindex(field, i, j, k) returns field.data[i, j, k]
broadcast follows this pattern:

julia> Meta.@lower c1[1, :, 1] .= c2[2, :, 3]
:($(Expr(:thunk, CodeInfo(
    @ none within `top-level scope`
1 ─ %1 = Base.dotview(c1, 1, :, 1)
│   %2 = Base.getindex(c2, 2, :, 3)
│   %3 = Base.broadcasted(Base.identity, %2)
│   %4 = Base.materialize!(%1, %3)
└──      return %4
))))

and

julia> axes(Base.dotview(c1, 1, :, 1))
(1:1, Base.OneTo(3), 1:1)

julia> axes(Base.getindex(c2, 2, :, 3))
(OffsetArrays.IdOffsetRange(values=-2:6, indices=-2:6),)

This will work:

c1[1, :, 1] .= view(c2, 1, :, 1)

In general c1[1, :, 1] .= c2[1, :, 1] is not really GPU-compatible. Anyways, if we want to support this behavior we need to change getindex(field, i, j, k) to return a view instead of the underlying data in case of a Colon index.
We might need to think this a bit though, because it would be a big change and might break other implementations

simone-silvestri · 2023-10-25T01:36:42Z

This in src/Fields/field.jl should allow that type of broadcasting

@propagate_inbounds Base.getindex(f::Field, ::Colon, j, k) = view(f, :, j, k)
@propagate_inbounds Base.getindex(f::Field, i, ::Colon, k) = view(f, i, :, k)
@propagate_inbounds Base.getindex(f::Field, i, j, ::Colon) = view(f, i, j, :)
@propagate_inbounds Base.getindex(f::Field, ::Colon, ::Colon, k) = view(f, :, :, k)
@propagate_inbounds Base.getindex(f::Field, ::Colon, j, ::Colon) = view(f, :, j, :)
@propagate_inbounds Base.getindex(f::Field, i, ::Colon, ::Colon) = view(f, i, :, :)
@propagate_inbounds Base.getindex(f::Field, ::Colon, ::Colon, ::Colon) = view(f, :, :, :)

navidcy · 2023-10-25T02:02:24Z

I see!

is this

c1[1, :, 1] .= view(c2, 1, :, 1)

or this

view(c1, 1, :, 1) .= view(c2, 1, :, 1)

preferred?

navidcy · 2023-10-25T02:04:03Z

This in src/Fields/field.jl should allow that type of broadcasting

@propagate_inbounds Base.getindex(f::Field, ::Colon, j, k) = view(f, :, j, k)
@propagate_inbounds Base.getindex(f::Field, i, ::Colon, k) = view(f, i, :, k)
@propagate_inbounds Base.getindex(f::Field, i, j, ::Colon) = view(f, i, j, :)
@propagate_inbounds Base.getindex(f::Field, ::Colon, ::Colon, k) = view(f, :, :, k)
@propagate_inbounds Base.getindex(f::Field, ::Colon, j, ::Colon) = view(f, :, j, :)
@propagate_inbounds Base.getindex(f::Field, i, ::Colon, ::Colon) = view(f, i, :, :)
@propagate_inbounds Base.getindex(f::Field, ::Colon, ::Colon, ::Colon) = view(f, :, :, :)

Would this be a good idea to add?

simone-silvestri · 2023-10-25T02:08:43Z

I see!

is this

c1[1, :, 1] .= view(c2, 1, :, 1)

or this

view(c1, 1, :, 1) .= view(c2, 1, :, 1)

preferred?

The second one is preferred because it's GPU-compatible.
I am not sure if we should add that behavior. That would essentially change the underlying functionality of getindex to return different types when using different indices, so it would be a big change. I am not sure allowing sliced broadcasts is enough to justify this change.

simone-silvestri · 2023-10-25T02:09:29Z

Actually, I have to rectify, both commands do the same thing because c1[1, :, 1] on the LHS calls view(c1, 1, :, 1)

navidcy · 2023-10-25T02:10:46Z

Ok cool. Gotcha. I’m closing this.

glwagner · 2023-10-25T02:17:16Z

I see!
is this
c1[1, :, 1] .= view(c2, 1, :, 1)
or this
view(c1, 1, :, 1) .= view(c2, 1, :, 1)
preferred?
The second one is preferred because it's GPU-compatible. I am not sure if we should add that behavior. That would essentially change the underlying functionality of getindex to return different types when using different indices, so it would be a big change. I am not sure allowing sliced broadcasts is enough to justify this change.

I think those are actually identical (or they could be --- for arrays broadcasting will call view under the hood)

siddharthabishnu · 2024-04-02T00:11:50Z

The issue reappears when view is used with reverse. Consider the MWE:

julia> using Oceananigans
julia> grid = ConformalCubedSphereGrid(; z=(-1, 0), panel_size=(4, 4, 1), horizontal_direction_halo = 4, z_halo = 1)
julia> u = Field{Face, Center, Center}(grid)
julia> set!(u, 1)
julia> region, region_E, region_W, Nc, Hc, k, plmn = 1, 2, 5, grid.Nx, grid.Hx, 1, -1

Now, typing

julia> u[region][1, Nc+1:Nc+Hc, k] = reverse(u[region_W][1, Nc+1-Hc:Nc, k]) * plmn

works, but typing

julia> u[region][1, Nc+1:Nc+Hc, k] .= reverse(view(u[region_W], 1, Nc+1-Hc:Nc, k)) * plmn

throws the error

ERROR: MethodError: no method matching isinteger(::CartesianIndices{3, Tuple{UnitRange{Int64}, UnitRange{Int64}, UnitRange{Int64}}})

Closest candidates are:
  isinteger(::BigFloat)
   @ Base mpfr.jl:960
  isinteger(::Missing)
   @ Base missing.jl:101
  isinteger(::FixedPointNumbers.FixedPoint)
   @ FixedPointNumbers ~/.julia/packages/FixedPointNumbers/HAGk2/src/FixedPointNumbers.jl:101
  ...

Stacktrace:
  [1] validate_index(idx::CartesianIndices{3, Tuple{UnitRange{Int64}, UnitRange{Int64}, UnitRange{Int64}}}, loc::Face, topo::FullyConnected, N::Int64, H::Int64)
    @ Oceananigans.Grids /Users/Sid/Library/CloudStorage/Dropbox/StudyFolder/PostDocMITDesktop/Codes/Oceananigans/cs-grid-metrics/src/Grids/input_validation.jl:196
  [2] map
    @ ./tuple.jl:342 [inlined]
  [3] validate_indices(indices::Tuple{…}, loc::Tuple{…}, topo::Tuple{…}, sz::Tuple{…}, halo_sz::Tuple{…})
    @ Oceananigans.Grids /Users/Sid/Library/CloudStorage/Dropbox/StudyFolder/PostDocMITDesktop/Codes/Oceananigans/cs-grid-metrics/src/Grids/input_validation.jl:214
  [4] validate_indices(indices::Tuple{…}, loc::Tuple{…}, grid::ZRegOrthogonalSphericalShellGrid{…})
    @ Oceananigans.Grids /Users/Sid/Library/CloudStorage/Dropbox/StudyFolder/PostDocMITDesktop/Codes/Oceananigans/cs-grid-metrics/src/Grids/input_validation.jl:211
  [5] view(f::Field{…}, i::CartesianIndices{…}, j::Function, k::Function)
    @ Oceananigans.Fields /Users/Sid/Library/CloudStorage/Dropbox/StudyFolder/PostDocMITDesktop/Codes/Oceananigans/cs-grid-metrics/src/Fields/field.jl:308
  [6] view(f::Field{…}, i::CartesianIndices{…})
    @ Oceananigans.Fields /Users/Sid/Library/CloudStorage/Dropbox/StudyFolder/PostDocMITDesktop/Codes/Oceananigans/cs-grid-metrics/src/Fields/field.jl:339
  [7] _reverse!(A::Field{…}, dims::Tuple{…})
    @ Base ./arraymath.jl:95
  [8] _reverse!
    @ ./arraymath.jl:71 [inlined]
  [9] #reverse!#274
    @ ./arraymath.jl:70 [inlined]
 [10] _reverse(A::Field{…}, dims::Function)
    @ Base ./arraymath.jl:60
 [11] reverse(A::Field{Face, Center, Center, Nothing, ZRegOrthogonalSphericalShellGrid{…}, Tuple{…}, OffsetArray{…}, Float64, FieldBoundaryConditions{…}, Nothing, Oceananigans.Fields.FieldBoundaryBuffers{…}})
    @ Base ./arraymath.jl:59
 [12] top-level scope
    @ REPL[199]:1
Some type information was truncated. Use `show(err)` to see complete types.

If the vector is formed by extracting multiple elements from the first dimension of the field on the RHS (as opposed to the second dimension as above), e.g.,

julia> u[region][Nc+1, 1-Hc:0, k] .= reverse(view(u[region_E], 2:Hc+1, 1, k))

no error message pops up but u[region][Nc+1, 1-Hc:0, k] is filled with junk values as shown below:

1×4×1 Field{Face, Center, Center} on OrthogonalSphericalShellGrid on CPU
├── grid: 4×4×1 OrthogonalSphericalShellGrid{Float64, FullyConnected, FullyConnected, Bounded} on CPU with 4×4×1 halo and with precomputed metrics
├── boundary conditions: FieldBoundaryConditions
│   └── west: Nothing, east: Nothing, south: Nothing, north: Nothing, bottom: Nothing, top: Nothing, immersed: ZeroFlux
├── indices: (5:5, -3:0, 1:1)
└── data: 1×4×1 OffsetArray(view(::Array{Float64, 3}, 9:9, 1:4, 2:2), 5:5, -3:0, 1:1) with eltype Float64 with indices 5:5×-3:0×1:1
    └── max=2.1234e-314, min=5.0e-324, mean=5.30853e-315

glwagner · 2024-04-02T05:59:51Z

This in src/Fields/field.jl should allow that type of broadcasting

@propagate_inbounds Base.getindex(f::Field, ::Colon, j, k) = view(f, :, j, k)
@propagate_inbounds Base.getindex(f::Field, i, ::Colon, k) = view(f, i, :, k)
@propagate_inbounds Base.getindex(f::Field, i, j, ::Colon) = view(f, i, j, :)
@propagate_inbounds Base.getindex(f::Field, ::Colon, ::Colon, k) = view(f, :, :, k)
@propagate_inbounds Base.getindex(f::Field, ::Colon, j, ::Colon) = view(f, :, j, :)
@propagate_inbounds Base.getindex(f::Field, i, ::Colon, ::Colon) = view(f, i, :, :)
@propagate_inbounds Base.getindex(f::Field, ::Colon, ::Colon, ::Colon) = view(f, :, :, :)

Would this be a good idea to add?

Don't do this

glwagner · 2024-04-02T06:10:53Z

Hmm. What is

view(u[region_W], 1, Nc+1-Hc:Nc, k)

is this a Field? Next question, what is

reverse(view(u[region_W], 1, Nc+1-Hc:Nc, k))

If the former is a Field I don't think we've defined reverse on Field. But you could do that...

siddharthabishnu · 2024-04-03T00:23:00Z

Hmm. What is
view(u[region_W], 1, Nc+1-Hc:Nc, k)
is this a Field? Next question, what is
reverse(view(u[region_W], 1, Nc+1-Hc:Nc, k))
If the former is a Field I don't think we've defined reverse on Field. But you could do that...

@glwagner, yes, view(u[region_W], 1, Nc+1-Hc:Nc, k) is a WindowedField, on which reverse is not defined yet. But reverse(view(u[region_W], 1, Nc+1-Hc:Nc, k).data) is a valid operation. However, simply typing

u[region][1, Nc+1:Nc+Hc, k] .= reverse(view(u[region_W], 1, Nc+1-Hc:Nc, k).data) * plmn
u[region][Nc+1, 1-Hc:0, k] .= reverse(view(u[region_E], 2:Hc+1, 1, k))

might throw the above-mentioned error due to a dimension mismatch or fill the LHS with junk or incorrect values. So, a safe solution is

view(u[region], 1, Nc+1:Nc+Hc, k).data .= reshape(reverse(view(u[region_W], 1, Nc+1-Hc:Nc, k).data) * plmn, 1:1, Nc+1:Nc+Hc, k:k)
view(u[region], Nc+1, 1-Hc:0, k).data .= reshape(reverse(view(u[region_E], 2:Hc+1, 1, k).data), Nc+1:Nc+1, 1-Hc:0, k:k)

@glwagner, @navidcy, if you guys are fine with this fix, I will close the issue.

navidcy · 2024-04-03T04:22:31Z

might either throw an error or fill the LHS with junk or incorrect values

I'm not that happy mostly because it seems we don't understand? Why "might"?

siddharthabishnu · 2024-04-03T13:33:20Z

might either throw an error or fill the LHS with junk or incorrect values

I'm not that happy mostly because it seems we don't understand? Why "might"?

I guess I used "might" because I don't fully understand it either. What I meant to say is that the line

u[region][1, Nc+1:Nc+Hc, k] .= reverse(view(u[region_W], 1, Nc+1-Hc:Nc, k).data) * plmn

throws an error due to a dimension mismatch and the line

u[region][Nc+1, 1-Hc:0, k] .= reverse(view(u[region_E], 2:Hc+1, 1, k))

fills the LHS with junk values. Hence the workaround I proposed above.

siddharthabishnu · 2024-04-03T13:57:24Z

Here's another observation:

u[region][1, 1-Hc:0, k] .= view(u[region_W], Nc+1-Hc:Nc, 1, k)

neither updates the LHS with any value nor triggers an error. The fix is either

u[region][1, 1-Hc:0, k] = u[region_W][Nc+1-Hc:Nc, 1, k]

or

view(u[region], 1, 1-Hc:0, k).data .= reshape(view(u[region_W], Nc+1-Hc:Nc, 1, k).data, 1:1, 1-Hc:0, k:k)

(preferred).

glwagner · 2024-04-03T13:59:56Z

@siddharthabishnu if anything that solution isn't quite satisfactory because it is hard to read.

Is there a reason why you can't use parent rather than manually writing .data? There's also a function data(::Field) that you can use if you need. I suspect though that in this case either using parent(f) or data(f) will yield the same result.

Can you print the REPL output of

view(u[region_W], Nc+1-Hc:Nc, 1, k)

I'm still curious what reverse(f::Field) outputs too.

siddharthabishnu · 2024-04-03T14:11:16Z

@siddharthabishnu if anything that solution isn't quite satisfactory because it is hard to read.

Is there a reason why you can't use parent rather than manually writing .data? There's also a function data(::Field) that you can use if you need. I suspect though that in this case either using parent(f) or data(f) will yield the same result.

Can you print the REPL output of
view(u[region_W], Nc+1-Hc:Nc, 1, k)
I'm still curious what reverse(f::Field) outputs too.

@glwagner, view(u[region_W], Nc+1-Hc:Nc, 1, k) results in a WindowedField. The REPL output is

julia> view(u[region_W], Nc+1-Hc:Nc, 1, k)
4×1×1 Field{Face, Center, Center} on OrthogonalSphericalShellGrid on CPU
├── grid: 4×4×1 OrthogonalSphericalShellGrid{Float64, FullyConnected, FullyConnected, Bounded} on CPU with 4×4×1 halo and with precomputed metrics
├── boundary conditions: FieldBoundaryConditions
│   └── west: Nothing, east: Nothing, south: Nothing, north: Nothing, bottom: Nothing, top: Nothing, immersed: ZeroFlux
├── indices: (1:4, 1:1, 1:1)
└── data: 4×1×1 OffsetArray(view(::Array{Float64, 3}, 5:8, 5:5, 2:2), 1:4, 1:1, 1:1) with eltype Float64 with indices 1:4×1:1×1:1
    └── max=1.0, min=1.0, mean=1.0

Also, reverse(u[region][1, 1:Nc, k]) outputs a vector, as shown below:

julia> reverse(u[region][1, 1:Nc, k])
4-element Vector{Float64}:
 1.0
 1.0
 1.0
 1.0

I experimented with the parent function a while back but cannot recall the specific issues encountered. I will revisit this and share my findings here.

navidcy · 2024-04-03T14:24:34Z

@glwagner is curious what reverse(f::Field) gives though; you showed them what reverse(a::Array) does.

navidcy · 2024-04-03T14:47:18Z

look..

julia> typeof(u[1])
Field{Face, Center, Center, Nothing, ZRegOrthogonalSphericalShellGrid{Float64, FullyConnected, FullyConnected, Bounded, OffsetMatrix{Float64, Matrix{Float64}}, OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, Float64, @NamedTuple{ξ::Tuple{Float64, Float64}, η::Tuple{Float64, Float64}, rotation::Rotations.RotXY{Float64}}, CPU}, Tuple{Colon, Colon, Colon}, OffsetArray{Float64, 3, Array{Float64, 3}}, Float64, FieldBoundaryConditions{BoundaryCondition{Oceananigans.BoundaryConditions.MultiRegionCommunication, Oceananigans.MultiRegion.CubedSphereRegionalConnectivity{West, North, Oceananigans.MultiRegion.↻}}, BoundaryCondition{Oceananigans.BoundaryConditions.MultiRegionCommunication, Oceananigans.MultiRegion.CubedSphereRegionalConnectivity{East, West, Nothing}}, BoundaryCondition{Oceananigans.BoundaryConditions.MultiRegionCommunication, Oceananigans.MultiRegion.CubedSphereRegionalConnectivity{South, North, Nothing}}, BoundaryCondition{Oceananigans.BoundaryConditions.MultiRegionCommunication, Oceananigans.MultiRegion.CubedSphereRegionalConnectivity{North, West, Oceananigans.MultiRegion.↺}}, BoundaryCondition{Flux, Nothing}, BoundaryCondition{Flux, Nothing}, BoundaryCondition{Flux, Nothing}}, Nothing, Oceananigans.Fields.FieldBoundaryBuffers{@NamedTuple{send::Array{Float64, 3}, recv::Array{Float64, 3}}, @NamedTuple{send::Array{Float64, 3}, recv::Array{Float64, 3}}, @NamedTuple{send::Array{Float64, 3}, recv::Array{Float64, 3}}, @NamedTuple{send::Array{Float64, 3}, recv::Array{Float64, 3}}, Nothing, Nothing, Nothing, Nothing}}

julia> @which reverse(u[1])
reverse(A::AbstractArray; dims)
     @ Base arraymath.jl:59

isn't it strange that it uses the reverse for AbstractArray?

glwagner · 2024-04-03T18:03:59Z

Hmm, Sid I think you made a mistake because

u[region][1, 1:Nc, k]

is not the same thing as view(u[region], 1, 1:Nc, k). Right? Earlier, you were trying to use reverse(view(u[region_W], 1, Nc+1-Hc:Nc, k)).

@navidcy Field is an AbstractArray, because

Oceananigans.jl/src/Fields/abstract_field.jl

Line 29 in f7042b7

 abstract type AbstractField{LX, LY, LZ, G <: GridOrNothing, T, N} <: AbstractArray{T, N} end 

So, reverse is defined. But it might not be working properly, let's look into it. If it doesn't work, we could define it. Or, we can use

slice = view(u[region], 1, 1:Nc, k)
reversed_slice = reverse(parent(slice))

glwagner · 2024-04-03T18:04:14Z

I think reverse allocates memory though, so I don't know if we want to use it...

navidcy · 2024-04-03T18:18:16Z

True , we should avoid reverse.

glwagner · 2024-04-03T18:37:24Z

Well, it looks like reverse does execute (on CPU, but probably not on GPU if I had to guess). But its not correct:

julia> grid = RectilinearGrid(size=(2, 3, 4), x=(0, 1), y=(0, 1), z=(0, 1));

julia> c = CenterField(grid);

julia> set!(c, (x, y, z) -> randn())
2×3×4 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 2×3×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
└── data: 8×9×10 OffsetArray(::Array{Float64, 3}, -2:5, -2:6, -2:7) with eltype Float64 with indices -2:5×-2:6×-2:7
    └── max=2.78883, min=-2.29288, mean=0.1837

julia> interior(c)
2×3×4 view(::Array{Float64, 3}, 4:5, 4:6, 4:7) with eltype Float64:
[:, :, 1] =
 0.597135  -0.145186  0.705053
 1.64813   -0.927427  1.26691

[:, :, 2] =
  0.477191  -0.34755   -1.52352
 -0.238858   0.556976  -0.554776

[:, :, 3] =
 -2.29288   1.09366   -1.42584
  2.32618  -0.692534   1.32913

[:, :, 4] =
  0.411535  -0.418008  0.754811
 -0.52217   -0.458004  2.78883

julia> rc = reverse(c);

julia> interior(rc)
2×3×4 view(::Array{Float64, 3}, 4:5, 4:6, 4:7) with eltype Float64:
[:, :, 1] =
 0.0       0.0       0.0
 1.64813  -0.927427  1.26691

[:, :, 2] =
  0.0       0.0        0.0
 -0.238858  0.556976  -0.554776

[:, :, 3] =
 0.0       0.0       0.0
 2.32618  -0.692534  1.32913

[:, :, 4] =
  0.0       0.0       0.0
 -0.52217  -0.458004  2.78883

glwagner · 2024-04-03T18:39:03Z

For matrices, I think this reverses the matrix as if it were unrolled by linear indexing. So:

julia> a = rand(2, 3)
2×3 Matrix{Float64}:
 0.0277495  0.500604  0.304882
 0.146228   0.71934   0.705938

julia> a[:]
6-element Vector{Float64}:
 0.027749521297110946
 0.14622819860123326
 0.5006038914167457
 0.7193398558177163
 0.30488229025160674
 0.7059383732478947

julia> reverse(a)[:]
6-element Vector{Float64}:
 0.7059383732478947
 0.30488229025160674
 0.7193398558177163
 0.5006038914167457
 0.14622819860123326
 0.027749521297110946

glwagner · 2024-04-03T18:47:49Z

Probably not the best use of time to try to get reverse to work. But if someone is interested in looking further, we can see that reverse starts by building a new object with copymutable:

https://github.com/JuliaLang/julia/blob/d7dc9a8cc8f2aebf04d5cecc8625be250169644b/base/arraymath.jl#L59-L60

copymutable works correctly because it relies on 1) similar and 2) copyto!, both of which we define for Field. So that's nice, showing the power of extending those functions.

The error that @siddharthabishnu finds occurs because we haven't defined view for CartesianIndices. This is used by _reverse:

https://github.com/JuliaLang/julia/blob/d7dc9a8cc8f2aebf04d5cecc8625be250169644b/base/arraymath.jl#L95

The fix is pretty easy I think. We simply need to extend reverse!:

Base.reverse!(f::Field; dims=:) = reverse!(parent(f); dims)

This seems to work as expected:

julia> reverse!(parent(c)); interior(c)[:]
6-element Vector{Float64}:
 -0.14945052605940137
 -0.25420686637327716
  0.36399517922432356
 -0.03994475992128757
  0.018968176930285508
 -0.5607261175546058

julia> reverse!(parent(c)); interior(c)[:]
6-element Vector{Float64}:
 -0.5607261175546058
  0.018968176930285508
 -0.03994475992128757
  0.36399517922432356
 -0.25420686637327716
 -0.14945052605940137

But again I don't know why we want reverse for Field. Do we?

glwagner · 2024-04-03T18:52:02Z

I'm sort of amazed but it looks like this can work for windowed fields too:

julia> cv = view(c, 1, :, 2:3)

julia> reverse!(parent(cv)); interior(cv)[:]
4-element Vector{Float64}:
 -0.03994475992128757
  0.36399517922432356
 -0.25420686637327716
 -0.14945052605940137

julia> reverse!(parent(cv)); interior(cv)[:]
4-element Vector{Float64}:
 -0.14945052605940137
 -0.25420686637327716
  0.36399517922432356
 -0.03994475992128757

So apparently reverse! works properly for SubArray. And that's what parent(cv) is:

julia> parent(cv)
1×4×2 view(::Array{Float64, 3}, 2:2, :, 3:4) with eltype Float64:
[:, :, 1] =
 0.0  -0.0399448  0.363995  0.0

[:, :, 2] =
 0.0  -0.254207  -0.149451  0.0

Doubtful this all works on the GPU though, but who knows.

navidcy · 2024-04-03T18:54:51Z

We don’t need reverse for Fields. We need to go away from reverse since it allocates

glwagner · 2024-04-03T18:57:39Z

Also this isn't related to broadcasting is it?

siddharthabishnu · 2024-04-04T00:51:17Z

@glwagner and @navidcy, thanks for looking into this issue, and providing the necessary clarifications and insights into the reverse operator. It's really helpful.

@glwagner, let's now take a step back from reverse. I agree u[1][1, 1:Nc, k] is not the same thing as view(u[1], 1, 1:Nc, k). u[1][1, 1:Nc, k] is a vector and view(u[1], 1, 1:Nc, k) is a WindowedField.

julia> u[1][1, 1:Nc, k]
4-element Vector{Float64}:
 1.0
 1.0
 1.0
 1.0

julia> view(u[1], 1, 1:Nc, k)
1×4×1 Field{Face, Center, Center} on OrthogonalSphericalShellGrid on CPU
├── grid: 4×4×1 OrthogonalSphericalShellGrid{Float64, FullyConnected, FullyConnected, Bounded} on CPU with 4×4×1 halo and with precomputed metrics
├── boundary conditions: FieldBoundaryConditions
│   └── west: Nothing, east: Nothing, south: Nothing, north: Nothing, bottom: Nothing, top: Nothing, immersed: ZeroFlux
├── indices: (1:1, 1:4, 1:1)
└── data: 1×4×1 OffsetArray(view(::Array{Float64, 3}, 5:5, 5:8, 2:2), 1:1, 1:4, 1:1) with eltype Float64 with indices 1:1×1:4×1:1
    └── max=1.0, min=1.0, mean=1.0

However, the outcome of

u[1][1:Nc, 1, k] .= u[5][2:Nc+1, 2, k]
u[1][1:Nc, 1, k] .= view(u[5], 2:Nc+1, 2, k)
view(u[1], 1:Nc, 1, k) .= u[5][2:Nc+1, 2, k]
view(u[1], 1:Nc, 1, k) .= view(u[5], 2:Nc+1, 2, k)
is the same as shown below:

julia> set!(u[1], 1)
4×4×1 Field{Face, Center, Center} on OrthogonalSphericalShellGrid on CPU
├── grid: 4×4×1 OrthogonalSphericalShellGrid{Float64, FullyConnected, FullyConnected, Bounded} on CPU with 4×4×1 halo and with precomputed metrics
├── boundary conditions: FieldBoundaryConditions
│   └── west: MultiRegionCommunication, east: MultiRegionCommunication, south: MultiRegionCommunication, north: MultiRegionCommunication, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
└── data: 12×12×3 OffsetArray(::Array{Float64, 3}, -3:8, -3:8, 0:2) with eltype Float64 with indices -3:8×-3:8×0:2
    └── max=1.0, min=1.0, mean=1.0

julia> set!(u[5], 2)
4×4×1 Field{Face, Center, Center} on OrthogonalSphericalShellGrid on CPU
├── grid: 4×4×1 OrthogonalSphericalShellGrid{Float64, FullyConnected, FullyConnected, Bounded} on CPU with 4×4×1 halo and with precomputed metrics
├── boundary conditions: FieldBoundaryConditions
│   └── west: MultiRegionCommunication, east: MultiRegionCommunication, south: MultiRegionCommunication, north: MultiRegionCommunication, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
└── data: 12×12×3 OffsetArray(::Array{Float64, 3}, -3:8, -3:8, 0:2) with eltype Float64 with indices -3:8×-3:8×0:2
    └── max=2.0, min=2.0, mean=2.0

julia> u[1][1:Nc, 1, k] .= u[5][2:Nc+1, 2, k]
4×1×1 view(::Array{Float64, 3}, 5:8, 5:5, 2:2) with eltype Float64:
[:, :, 1] =
 2.0
 2.0
 2.0
 0.0

julia> u[1][:, :, k]
12×12 OffsetArray(::Matrix{Float64}, -3:8, -3:8) with eltype Float64 with indices -3:8×-3:8:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  2.0  1.0  1.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  2.0  1.0  1.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  2.0  1.0  1.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  1.0  1.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0

julia> set!(u[1], 1)
4×4×1 Field{Face, Center, Center} on OrthogonalSphericalShellGrid on CPU
├── grid: 4×4×1 OrthogonalSphericalShellGrid{Float64, FullyConnected, FullyConnected, Bounded} on CPU with 4×4×1 halo and with precomputed metrics
├── boundary conditions: FieldBoundaryConditions
│   └── west: MultiRegionCommunication, east: MultiRegionCommunication, south: MultiRegionCommunication, north: MultiRegionCommunication, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
└── data: 12×12×3 OffsetArray(::Array{Float64, 3}, -3:8, -3:8, 0:2) with eltype Float64 with indices -3:8×-3:8×0:2
    └── max=1.0, min=1.0, mean=1.0

julia> u[1][1:Nc, 1, k] .= view(u[5], 2:Nc+1, 2, k)
4×1×1 Field{Face, Center, Center} on OrthogonalSphericalShellGrid on CPU
├── grid: 4×4×1 OrthogonalSphericalShellGrid{Float64, FullyConnected, FullyConnected, Bounded} on CPU with 4×4×1 halo and with precomputed metrics
├── boundary conditions: FieldBoundaryConditions
│   └── west: Nothing, east: Nothing, south: Nothing, north: Nothing, bottom: Nothing, top: Nothing, immersed: ZeroFlux
├── indices: (1:4, 1:1, 1:1)
└── data: 4×1×1 OffsetArray(view(::Array{Float64, 3}, 5:8, 5:5, 2:2), 1:4, 1:1, 1:1) with eltype Float64 with indices 1:4×1:1×1:1
    └── max=2.0, min=2.0, mean=2.0

julia> u[1][:, :, k]
12×12 OffsetArray(::Matrix{Float64}, -3:8, -3:8) with eltype Float64 with indices -3:8×-3:8:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  2.0  1.0  1.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  2.0  1.0  1.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  2.0  1.0  1.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  2.0  1.0  1.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0

julia> set!(u[1], 1)
4×4×1 Field{Face, Center, Center} on OrthogonalSphericalShellGrid on CPU
├── grid: 4×4×1 OrthogonalSphericalShellGrid{Float64, FullyConnected, FullyConnected, Bounded} on CPU with 4×4×1 halo and with precomputed metrics
├── boundary conditions: FieldBoundaryConditions
│   └── west: MultiRegionCommunication, east: MultiRegionCommunication, south: MultiRegionCommunication, north: MultiRegionCommunication, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
└── data: 12×12×3 OffsetArray(::Array{Float64, 3}, -3:8, -3:8, 0:2) with eltype Float64 with indices -3:8×-3:8×0:2
    └── max=1.0, min=1.0, mean=1.0

julia> view(u[1], 1:Nc, 1, k) .= u[5][2:Nc+1, 2, k]
4×1×1 view(::Array{Float64, 3}, 5:8, 5:5, 2:2) with eltype Float64:
[:, :, 1] =
 2.0
 2.0
 2.0
 0.0

julia> u[1][:, :, k]
12×12 OffsetArray(::Matrix{Float64}, -3:8, -3:8) with eltype Float64 with indices -3:8×-3:8:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  2.0  1.0  1.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  2.0  1.0  1.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  2.0  1.0  1.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  1.0  1.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0

julia> set!(u[1], 1)
4×4×1 Field{Face, Center, Center} on OrthogonalSphericalShellGrid on CPU
├── grid: 4×4×1 OrthogonalSphericalShellGrid{Float64, FullyConnected, FullyConnected, Bounded} on CPU with 4×4×1 halo and with precomputed metrics
├── boundary conditions: FieldBoundaryConditions
│   └── west: MultiRegionCommunication, east: MultiRegionCommunication, south: MultiRegionCommunication, north: MultiRegionCommunication, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
└── data: 12×12×3 OffsetArray(::Array{Float64, 3}, -3:8, -3:8, 0:2) with eltype Float64 with indices -3:8×-3:8×0:2
    └── max=1.0, min=1.0, mean=1.0

julia> view(u[1], 1:Nc, 1, k) .= view(u[5], 2:Nc+1, 2, k)
4×1×1 Field{Face, Center, Center} on OrthogonalSphericalShellGrid on CPU
├── grid: 4×4×1 OrthogonalSphericalShellGrid{Float64, FullyConnected, FullyConnected, Bounded} on CPU with 4×4×1 halo and with precomputed metrics
├── boundary conditions: FieldBoundaryConditions
│   └── west: Nothing, east: Nothing, south: Nothing, north: Nothing, bottom: Nothing, top: Nothing, immersed: ZeroFlux
├── indices: (1:4, 1:1, 1:1)
└── data: 4×1×1 OffsetArray(view(::Array{Float64, 3}, 5:8, 5:5, 2:2), 1:4, 1:1, 1:1) with eltype Float64 with indices 1:4×1:1×1:1
    └── max=2.0, min=2.0, mean=2.0

julia> u[1][:, :, k]
12×12 OffsetArray(::Matrix{Float64}, -3:8, -3:8) with eltype Float64 with indices -3:8×-3:8:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  2.0  1.0  1.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  2.0  1.0  1.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  2.0  1.0  1.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  2.0  1.0  1.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0

This is because in broadcast assignment operations, such as u[1][1:Nc, 1, k] .= u[5][2:Nc+1, 2, k], the data of the WindowedField is assigned. So, where's the issue? It's when we try any of the following operations:

u[1][1, 1:Nc, k] .= u[5][2:Nc+1, 2, k]
view(u[1], 1, 1:Nc, k) .= u[5][2:Nc+1, 2, k]
u[1][1, 1:Nc, k] .= u[5][2, 2:Nc+1, k]
view(u[1], 1, 1:Nc, k) .= u[5][2, 2:Nc+1, k]

Here the data extracted from the WindowedField on the RHS is always in the form of a column vector resulting in a dimension mismatch error. As a potential solution, we can consider assigning a view of the WindowedField as

u[1][1, 1:Nc, k] .= view(u[5], 2:Nc+1, 2, k)
view(u[1], 1, 1:Nc, k) .= view(u[5], 2:Nc+1, 2, k)
u[1][1, 1:Nc, k] .= view(u[5], 2, 2:Nc+1, k)
view(u[1], 1, 1:Nc, k) .= view(u[5], 2, 2:Nc+1, k)

or the transpose of the WindowedField as

u[1][1, 1:Nc, k] .= u[5][2:Nc+1, 2, k]'
view(u[1], 1, 1:Nc, k) .= u[5][2:Nc+1, 2, k]'
u[1][1, 1:Nc, k] .= u[5][2, 2:Nc+1, k]'
view(u[1], 1, 1:Nc, k) .= u[5][2, 2:Nc+1, k]'

It turns out only the second approach results in the correct solution

julia> u[1][:, :, 1]
12×12 OffsetArray(::Matrix{Float64}, -3:8, -3:8) with eltype Float64 with indices -3:8×-3:8:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  2.0  2.0  2.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  1.0  1.0  1.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  1.0  1.0  1.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  1.0  1.0  1.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0

whereas the first approach results in an incorrect solution

julia> u[1][:, :, 1]
12×12 OffsetArray(::Matrix{Float64}, -3:8, -3:8) with eltype Float64 with indices -3:8×-3:8:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  2.0  2.0  2.0  2.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  1.0  1.0  1.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  1.0  1.0  1.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  1.0  1.0  1.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0

I am not sure why that is the case. Even though view(u[5], 2, 2:Nc+1, k).data, view(u[5], 2, 2:Nc+1, k).data, view(u[5], 2:Nc+1, 2, k).data, or view(u[5], 2:Nc+1, 2, k).data results in the vector [2.0, 2.0, 2.0, 0.0] in row or column format, what is assigned is the vector [2.0, 2.0, 2.0, 2.0].

I believe this to be the root cause of the dimension mismatch bug and not the reverse operator. So, when reverse is applied to the transposed WindowedField (as in the second approach) when necessary, the correct solution is obtained without any issue.

glwagner · 2024-04-04T20:27:34Z

I think you're right that view(::Field) doesn't recover the same behavior as view(::Array). This is because Field are fixed to three dimensions; we don't have a way to express the concept of a Field with one dimension. So for example:

julia> using Oceananigans

julia> grid = RectilinearGrid(size=(2, 3, 4), x=(0, 1), y=(0, 1), z=(0, 1));

julia> c = CenterField(grid);

julia> size(view(c, 1, :, 2))
(1, 3, 1)

julia> size(view(parent(c), 1, :, 2))
(9,)

If our objective purely regards broadcasting, then we could probably generalize broadcasting so that Fields with singleton indices and/or Nothing location act more like reduced-dimensionality arrays. Is that what you would like to see?

glwagner · 2024-04-04T20:35:11Z

Rather than working in the complicated setting of a multiregion grid, it might be better to have a simpler context to look at the issue. Here's an MWE:

julia> grid = RectilinearGrid(size=(3, 3, 3), x=(0, 1), y=(0, 1), z=(0, 1));

julia> c = CenterField(grid); set!(c, (x, y, z) -> randn())
3×3×3 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 3×3×3 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
└── data: 9×9×9 OffsetArray(::Array{Float64, 3}, -2:6, -2:6, -2:6) with eltype Float64 with indices -2:6×-2:6×-2:6
    └── max=2.47269, min=-2.49748, mean=-0.0773234

julia> a = rand(3, 3, 3); b = rand(3, 3, 3);

Now we can illustrate how Field does not behave like Array. This works:

julia> view(a, :, 1, 1) .= view(b, 1, :, 1)
3-element view(::Array{Float64, 3}, :, 1, 1) with eltype Float64:
 0.23703924592859704
 0.02820508483630202
 0.7546236569353038

but this does not:

julia> view(a, :, 1, 1) .= view(c, 1, :, 1)

and instead throws DimensionMismatch. Does that boil down the issue @siddharthabishnu ?

Actually, it is curious that for some reason things do work out correctly if setting to a Field:

julia> d = CenterField(grid);

julia> view(d, :, 1, 1) .= view(c, 1, :, 1);

julia> interior(d, :, 1, 1)
3-element view(::Array{Float64, 3}, 4:6, 4, 4) with eltype Float64:
  0.06779701769558852
 -1.1932513091536738
 -1.1239966299690682

julia> interior(c, 1, :, 1)
3-element view(::Array{Float64, 3}, 4, 4:6, 4) with eltype Float64:
  0.06779701769558852
 -0.5935171905074843
  1.399692772040493

No idea why that works...

siddharthabishnu · 2024-04-05T01:35:51Z

Actually, it is curious that for some reason things do work out correctly if setting to a Field:

julia> d = CenterField(grid);

julia> view(d, :, 1, 1) .= view(c, 1, :, 1);

julia> interior(d, :, 1, 1)
3-element view(::Array{Float64, 3}, 4:6, 4, 4) with eltype Float64:
  0.06779701769558852
 -1.1932513091536738
 -1.1239966299690682

julia> interior(c, 1, :, 1)
3-element view(::Array{Float64, 3}, 4, 4:6, 4) with eltype Float64:
  0.06779701769558852
 -0.5935171905074843
  1.399692772040493

No idea why that works...

@glwagner, Thanks for creating the MWE without the complication of a multiregion grid and boiling down the issue. However, I disagree on your last point. Things do not work correctly in case of a Field. Yes, it doesn't throw an error, but the data of interior(c, 1, :, 1) is not assigned to interior(d, :, 1, 1) except for the first element (as seen in your MWE and mine below).

julia> a = rand(3, 3, 3); b = rand(3, 3, 3);

julia> grid = RectilinearGrid(size=(3, 3, 3), x=(0, 1), y=(0, 1), z=(0, 1));

julia> c = CenterField(grid); d = CenterField(grid);

julia> set!(c,a); set!(d,b);

julia> interior(d, :, 1, 1)
3-element view(::Array{Float64, 3}, 4:6, 4, 4) with eltype Float64:
 0.006219191907858468
 0.8963174540626572
 0.06824259256136866

julia> interior(c, 1, :, 1)
3-element view(::Array{Float64, 3}, 4, 4:6, 4) with eltype Float64:
 0.45223462463057673
 0.18046843834182114
 0.9782543808912043

julia> view(d, :, 1, 1) .= view(c, 1, :, 1)
3×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 3×3×3 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Nothing, north: Nothing, bottom: Nothing, top: Nothing, immersed: ZeroFlux
├── indices: (:, 1:1, 1:1)
└── data: 9×1×1 OffsetArray(view(::Array{Float64, 3}, :, 4:4, 4:4), -2:6, 1:1, 1:1) with eltype Float64 with indices -2:6×1:1×1:1
    └── max=0.697903, min=0.452235, mean=0.568357

julia> interior(d, :, 1, 1)
3-element view(::Array{Float64, 3}, 4:6, 4, 4) with eltype Float64:
 0.45223462463057673
 0.6979033041564312
 0.5549323977919884

However, the following works for the interior of the Field:

julia> interior(d, :, 1, 1) .= interior(c, 1, :, 1);

julia> interior(d, :, 1, 1)
3-element view(::Array{Float64, 3}, 4:6, 4, 4) with eltype Float64:
 0.45223462463057673
 0.18046843834182114
 0.9782543808912043

Moreover, since view(d, :, 1, 1).data is a 9x1x1 OffsetArray and view(c, 1, :, 1).data is a 1x9x1 OffsetArray,
view(d, :, 1, 1).data .= view(c, 1, :, 1).data leads to a dimension mismatch error. However, the following approach gives the correct result.

julia> set!(c,a); set!(d,b);

julia> view(d, :, 1, 1).data .= reshape(view(c, 1, :, 1).data, -2:6, 1:1, 1:1)
9×1×1 OffsetArray(view(::Array{Float64, 3}, :, 4:4, 4:4), -2:6, 1:1, 1:1) with eltype Float64 with indices -2:6×1:1×1:1:
[:, :, 1] =
 0.0
 0.0
 0.0
 0.45223462463057673
 0.18046843834182114
 0.9782543808912043
 0.0
 0.0
 0.0

siddharthabishnu · 2024-04-05T02:12:06Z

If our objective purely regards broadcasting, then we could probably generalize broadcasting so that Fields with singleton indices and/or Nothing location act more like reduced-dimensionality arrays. Is that what you would like to see?

Yes, it would be nice to have that feature.

glwagner · 2024-04-05T05:25:16Z

Yes, it doesn't throw an error, but the data of interior(c, 1, :, 1) is not assigned to interior(d, :, 1, 1) except for the first element (as seen in your MWE and mine below).

Heh, true. That's a relief actually...

glwagner · 2024-04-05T05:26:10Z

If our objective purely regards broadcasting, then we could probably generalize broadcasting so that Fields with singleton indices and/or Nothing location act more like reduced-dimensionality arrays. Is that what you would like to see?

Yes, it would be nice to have that feature.

Do you want to have a crack at it? I don't use broadcasting with fields so much personally. It's a bit slow for some reason, which we have never quite figured out. Convenient for some things but not to be relied on unless we can solve the performance issue.

navidcy added the question 💭 No such thing as a stupid question label Oct 24, 2023

navidcy closed this as completed Oct 25, 2023

navidcy reopened this Oct 25, 2023

navidcy closed this as completed Oct 25, 2023

navidcy mentioned this issue Mar 28, 2024

Fill halos of conformal cubed sphere grid coordinates and metrics in a single pass #3488

Merged

siddharthabishnu reopened this Apr 2, 2024

Broadcasting slices of fields doesn't work as expected #3361

Broadcasting slices of fields doesn't work as expected #3361

Comments

navidcy commented Oct 24, 2023 • edited

simone-silvestri commented Oct 25, 2023 • edited

simone-silvestri commented Oct 25, 2023

navidcy commented Oct 25, 2023 • edited

navidcy commented Oct 25, 2023

simone-silvestri commented Oct 25, 2023

simone-silvestri commented Oct 25, 2023

navidcy commented Oct 25, 2023

glwagner commented Oct 25, 2023

siddharthabishnu commented Apr 2, 2024 • edited

glwagner commented Apr 2, 2024

glwagner commented Apr 2, 2024

siddharthabishnu commented Apr 3, 2024 • edited

navidcy commented Apr 3, 2024

siddharthabishnu commented Apr 3, 2024 • edited

siddharthabishnu commented Apr 3, 2024 • edited

glwagner commented Apr 3, 2024

siddharthabishnu commented Apr 3, 2024

navidcy commented Apr 3, 2024

navidcy commented Apr 3, 2024

glwagner commented Apr 3, 2024

glwagner commented Apr 3, 2024

navidcy commented Apr 3, 2024

glwagner commented Apr 3, 2024

glwagner commented Apr 3, 2024

glwagner commented Apr 3, 2024

glwagner commented Apr 3, 2024

navidcy commented Apr 3, 2024

glwagner commented Apr 3, 2024

siddharthabishnu commented Apr 4, 2024 • edited

glwagner commented Apr 4, 2024 • edited

glwagner commented Apr 4, 2024 • edited

siddharthabishnu commented Apr 5, 2024 • edited

siddharthabishnu commented Apr 5, 2024

glwagner commented Apr 5, 2024

glwagner commented Apr 5, 2024

navidcy commented Oct 24, 2023 •

edited

simone-silvestri commented Oct 25, 2023 •

edited

navidcy commented Oct 25, 2023 •

edited

siddharthabishnu commented Apr 2, 2024 •

edited

siddharthabishnu commented Apr 3, 2024 •

edited

siddharthabishnu commented Apr 3, 2024 •

edited

siddharthabishnu commented Apr 3, 2024 •

edited

siddharthabishnu commented Apr 4, 2024 •

edited

glwagner commented Apr 4, 2024 •

edited

glwagner commented Apr 4, 2024 •

edited

siddharthabishnu commented Apr 5, 2024 •

edited