Bugfix for `parent_index_range` #3573

siddharthabishnu · 2024-05-01T04:32:24Z

Closes #3572.

src/Fields/field.jl

glwagner · 2024-05-01T17:29:47Z

There are a couple issues dedicated to index issues with the free surface that might be related. @siddharthabishnu can you try to track these down?

glwagner · 2024-05-01T17:31:46Z

Might be good to add a test (or more than just one) for windowed fields. Even a simple one like

windowed_field = CenterField(grid, indices=(:, :, 1:1))
@test view(windowed_field, :, :, 1:1) isa AbstractArray

or something.

@siddharthabishnu can you add that? Pretty basic test so probably belongs in test_field.jl:

https://github.com/CliMA/Oceananigans.jl/blob/sb/extend-parent-indices/test/test_field.jl

navidcy · 2024-05-02T04:49:09Z

A lot of tests seem to fail...

glwagner · 2024-05-02T15:00:11Z

A lot of tests seem to fail...

This is troubling

glwagner · 2024-05-02T15:03:59Z

Ok, I think it's because parent_index_range is used by the NetCDFOutputWriter. I suspect that it is being used incorrectly there though --- because I am nearly positive that this change does in fact fix a bug.

It's not the nicest code either

Oceananigans.jl/src/OutputWriters/netcdf_output_writer.jl

Lines 52 to 57 in ed77780

    
           Dict("xC" => parent(xnodes(grid, Center(); with_halos=true))[parent_index_range(indices["xC"][1], Center(), TX(), Hx)], 
        
                "xF" => parent(xnodes(grid, Face();   with_halos=true))[parent_index_range(indices["xF"][1],   Face(), TX(), Hx)], 
        
                "yC" => parent(ynodes(grid, Center(); with_halos=true))[parent_index_range(indices["yC"][2], Center(), TY(), Hy)], 
        
                "yF" => parent(ynodes(grid, Face();   with_halos=true))[parent_index_range(indices["yF"][2],   Face(), TY(), Hy)], 
        
                "zC" => parent(znodes(grid, Center(); with_halos=true))[parent_index_range(indices["zC"][3], Center(), TZ(), Hz)], 
        
                "zF" => parent(znodes(grid, Face();   with_halos=true))[parent_index_range(indices["zF"][3],   Face(), TZ(), Hz)])

glwagner · 2024-05-02T15:23:27Z

Ok I think I understand the issue. parent_index_range is being used for two different purposes.

In NetCDFOutputWriter the "parent" is an array that spans the entire dimension; we just want to translate indices defined on the grid to indices defined on the underlying parent array (I think the reason we need a function is because we need to treat Nothing locations specially).

In offset_windowed_data the indices argument are the indices of a windowed field itself.

I think we just need two different functions. I'll fix it.

glwagner · 2024-05-02T16:07:07Z

Ok here's a little more insight into the problem. The issue is actually with view, which doesn't work when called on windowed field. view can only create a windowed field from a full field, right now. In other words, view cannot act on a view. Haha.

I'm extending view to work correctly when the argument itself is already windowed.

This is probably the issue with outputting the free surface too. I can't remember exactly what went wrong there though. The output issue could also hint that things are more convoluted than they should be. Still though, we should be able to create views into fields that are already views. It should be fine.

navidcy · 2024-05-02T16:09:42Z

nice one!

glwagner · 2024-05-02T16:30:55Z

OY.

Ok I worked on it a bit. This works:

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

julia> c = CenterField(grid);

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

julia> cv = view(c, :, :, 2:3)
1×1×2 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 1×1×3 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: Nothing, top: Nothing, immersed: ZeroFlux
├── indices: (:, :, 2:3)
└── data: 3×3×2 OffsetArray(view(::Array{Float64, 3}, :, :, 5:6), 0:2, 0:2, 2:3) with eltype Float64 with indices 0:2×0:2×2:3
    └── max=0.92472, min=0.681611, mean=0.803165

julia> cvv = view(cv, :, :, 3)
1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 1×1×3 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: Nothing, top: Nothing, immersed: ZeroFlux
├── indices: (:, :, 3:3)
└── data: 3×3×1 OffsetArray(view(::Array{Float64, 3}, :, :, 6:6), 0:2, 0:2, 3:3) with eltype Float64 with indices 0:2×0:2×3:3
    └── max=0.681611, min=0.681611, mean=0.681611

julia> cvvv = view(cvv, :, :, 3)
1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 1×1×3 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: Nothing, top: Nothing, immersed: ZeroFlux
├── indices: (:, :, 3:3)
└── data: 3×3×1 OffsetArray(view(::Array{Float64, 3}, :, :, 6:6), 0:2, 0:2, 3:3) with eltype Float64 with indices 0:2×0:2×3:3
    └── max=0.681611, min=0.681611, mean=0.681611

julia> interior(c, :, :, 3) .= 0
1×1 view(::Array{Float64, 3}, 2:2, 2:2, 6) with eltype Float64:
 0.0

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

julia> cv
1×1×2 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 1×1×3 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: Nothing, top: Nothing, immersed: ZeroFlux
├── indices: (:, :, 2:3)
└── data: 3×3×2 OffsetArray(view(::Array{Float64, 3}, :, :, 5:6), 0:2, 0:2, 2:3) with eltype Float64 with indices 0:2×0:2×2:3
    └── max=0.92472, min=0.0, mean=0.46236

Note that one of the quirks of fields that come from views, is that parent returns a view:

julia> parent(cv)
3×3×2 view(::Array{Float64, 3}, :, :, 5:6) with eltype Float64:
[:, :, 1] =
 0.0  0.0      0.0
 0.0  0.92472  0.0
 0.0  0.0      0.0

[:, :, 2] =
 0.0  0.0  0.0
 0.0  0.0  0.0
 0.0  0.0  0.0

julia> parent(cvv)
3×3×1 view(::Array{Float64, 3}, :, :, 6:6) with eltype Float64:
[:, :, 1] =
 0.0  0.0  0.0
 0.0  0.0  0.0
 0.0  0.0  0.0

That's different than if you create a windowed field from scratch:

julia> b = CenterField(grid, indices=(:, :, 3))
1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 1×1×3 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: Nothing, top: Nothing, immersed: ZeroFlux
├── indices: (:, :, 3:3)
└── data: 3×3×1 OffsetArray(::Array{Float64, 3}, 0:2, 0:2, 3:3) with eltype Float64 with indices 0:2×0:2×3:3
    └── max=0.0, min=0.0, mean=0.0

julia> parent(b)
3×3×1 Array{Float64, 3}:
[:, :, 1] =
 0.0  0.0  0.0
 0.0  0.0  0.0
 0.0  0.0  0.0

This is how output works --- we use parent to get the data that needs to be saved down.

glwagner · 2024-05-02T16:35:20Z

@siddharthabishnu could you help by adding some simple unit tests for views of views? It'd be nice to test correctness, I think something like this would work:

grid = RectilinearGrid(size=(1, 1, 3), x=(0, 1), y=(0, 1), z=(0, 1))
c = CenterField(grid)
set!(c, (x, y, z) -> rand())

# First test that regular view is right
cv = view(c, :, :, 2:3)
@test c[1, 1, 2] == cv[1, 1, 2]
@test c[1, 1, 3] == cv[1, 1, 3]

# Now test views of views
cvv = view(cv, :, :, 3)
@test cv[1, 1, 3] == cvv[1, 1, 3]

We may also want to test that things error correctly like

julia> view(cv, :, :, 1)
ERROR: ArgumentError: view indices (Colon(), Colon(), 1) do not intersect field indices (Colon(), Colon(), 2:3)
Stacktrace:
 [1] view(f::Field{…}, i::Function, j::Function, k::Int64)
   @ Oceananigans.Fields ~/Projects/dev/Oceananigans.jl/src/Fields/field.jl:315
 [2] top-level scope
   @ REPL[25]:1
Some type information was truncated. Use `show(err)` to see complete types.

Maybe

@test_throws ArgumentError view(cv, :, :, 1)

Not sure what else we might want maybe @navidcy has ideas.

glwagner · 2024-05-02T16:36:00Z

Could do

@test_throws BoundsError cv[1, 1, 1]

glwagner · 2024-05-02T16:36:25Z

And to do those tests on the GPU we need CUDA.@allowscalar.

navidcy · 2024-05-02T16:53:53Z

Thanks for this @glwagner!

navidcy · 2024-05-02T17:03:22Z

I ran the doctests locally and this one:

Oceananigans.jl/src/Fields/field.jl

Lines 270 to 298 in f217b80

    
           ```jldoctest 
        
           julia> using Oceananigans 
        
           julia> grid = RectilinearGrid(size=(2, 3, 4), x=(0, 1), y=(0, 1), z=(0, 1)); 
        
           julia> c = 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 2×3×3 halo 
        
           ├── boundary conditions: FieldBoundaryConditions 
        
           │   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux 
        
           └── data: 6×9×10 OffsetArray(::Array{Float64, 3}, -1:4, -2:6, -2:7) with eltype Float64 with indices -1:4×-2:6×-2:7 
        
               └── max=0.0, min=0.0, mean=0.0 
        
           julia> c .= rand(size(c)...); 
        
           julia> v = view(c, :, 2:3, 1:2) 
        
           2×2×2 Field{Center, Center, Center} on RectilinearGrid on CPU 
        
           ├── grid: 2×3×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 2×3×3 halo 
        
           ├── boundary conditions: FieldBoundaryConditions 
        
           │   └── west: Periodic, east: Periodic, south: Nothing, north: Nothing, bottom: Nothing, top: Nothing, immersed: ZeroFlux 
        
           ├── indices: (:, 2:3, 1:2) 
        
           └── data: 6×2×2 OffsetArray(view(::Array{Float64, 3}, :, 5:6, 4:5), -1:4, 2:3, 1:2) with eltype Float64 with indices -1:4×2:3×1:2 
        
               └── max=0.972136, min=0.0149088, mean=0.59198 
        
           julia> size(v) 
        
           (2, 2, 2) 
        
           julia> v[2, 2, 2] == c[2, 2, 2] 
        
           true

fails.

I need to head out now but I wanted to document this.

navidcy · 2024-05-02T17:04:50Z

julia> using Oceananigans

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

julia> c = 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 2×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
└── data: 6×9×10 OffsetArray(::Array{Float64, 3}, -1:4, -2:6, -2:7) with eltype Float64 with indices -1:4×-2:6×-2:7
    └── max=0.0, min=0.0, mean=0.0

julia> c .= rand(size(c)...);

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

julia> size(v)
(2, 2, 2)

julia> v[2, 2, 2] == c[2, 2, 2]
false

seems like v is empty!

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

glwagner · 2024-05-02T18:02:44Z

Can't reproduce that @navidcy

julia> using Oceananigans
Precompiling Oceananigans
  1 dependency successfully precompiled in 11 seconds. 160 already precompiled.

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

julia> c = 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 2×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
└── data: 6×9×10 OffsetArray(::Array{Float64, 3}, -1:4, -2:6, -2:7) with eltype Float64 with indices -1:4×-2:6×-2:7
    └── max=0.0, min=0.0, mean=0.0

julia> c .= rand(size(c)...);

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

julia> size(v)
(2, 2, 2)

julia> v[2, 2, 2] == c[2, 2, 2]
true

navidcy · 2024-05-02T18:22:47Z

And all tests pass! Great! Ignore my messages!

glwagner · 2024-05-04T03:22:17Z

test/test_field.jl

+            @test all(cv[1, 1, i] == c[1, 1, i] for i in firstindex(interior(c), 3)+1:lastindex(interior(c), 3)-1)
+
+            # Now test the views of views
+            cvv = view(cv, :, :, firstindex(interior(c), 3)+2:lastindex(interior(c), 3)-2)


Maybe define k_tops = [Nz, Nz+1] and include that in the loop above

test/test_field.jl

glwagner · 2024-05-04T03:25:30Z

The tests use a lot of scalar indexing that's why they fail on the GPU:

https://buildkite.com/clima/oceananigans/builds/15604#018f40ed-787d-4e74-a5f4-ae1656fa3043/18-724

I think if we are comparing single numbers it makes sense to use @allowscalar.

If we are comparing vectors it could be nice to figure out how to get the tests to run without @allowscalar since presumably this is possible.

navidcy · 2024-05-08T06:27:01Z

@siddharthabishnu could you deal with this?

siddharthabishnu · 2024-05-09T20:41:34Z

The tests use a lot of scalar indexing that's why they fail on the GPU:

https://buildkite.com/clima/oceananigans/builds/15604#018f40ed-787d-4e74-a5f4-ae1656fa3043/18-724

I think if we are comparing single numbers it makes sense to use @allowscalar.

If we are comparing vectors it could be nice to figure out how to get the tests to run without @allowscalar since presumably this is possible.

Given that we're comparing elements of vectors with a maximum length of 6, I opted to use CUDA.@allowscalar. The impact on performance in this situation is minimal when running on a GPU.

siddharthabishnu · 2024-05-09T21:21:32Z

@siddharthabishnu could you deal with this?

Done. I think it's ready to merge once the tests pass.

siddharthabishnu · 2024-05-13T19:21:07Z

@glwagner, @navidcy, could one of you please review the changes and approve them if you are satisfied?

glwagner · 2024-05-13T20:40:43Z

test/test_field.jl

+                # @test_throws BoundsError cvv[:, :, 7:8] for a ZFaceField with bounded topology in z else @test_throws BoundsError cvv[:, :, 6:7].
+                @test_throws BoundsError cvvv[:, :, 1:1+2]
+                @test_throws BoundsError cvvv[:, :, k_top-2:k_top]
+                # @test_throws BoundsError cvvv[:, :, 6:8] for a ZFaceField with bounded topology in z else @test_throws BoundsError cvvv[:, :, 5:7].


@siddharthabishnu can you put @allowscalar around only the specific lines where it is needed? I thought by using k_top we eliminated the need for the comments. Not sure what the comments mean

@glwagner, agreed. Resolved in commits 7d68c02 and 427cb51.

glwagner · 2024-05-13T20:41:14Z

test/test_field.jl

+                @test_throws BoundsError cvv[:, :, k_top-1:k_top]
+                # @test_throws BoundsError cvv[:, :, 7:8] for a ZFaceField with bounded topology in z else @test_throws BoundsError cvv[:, :, 6:7].
+                @test_throws BoundsError cvvv[:, :, 1:1+2]
+                @test_throws BoundsError cvvv[:, :, k_top-2:k_top]


should this use interior? We don't ever directly index fields in the user interface

Suggested change

@test_throws BoundsError cvvv[:, :, k_top-2:k_top]

@test_throws BoundsError cvvv[:, :, k_top-2:k_top]

Unfortunately, not, unless we also extend the interior index range like we did with the parent index range. Consider the following MWE:

julia> using Oceananigans julia> Nx, Ny, Nz = 1, 1, 7; julia> grid = RectilinearGrid(size=(Nx, Ny, Nz), x=(0, 1), y=(0, 1), z=(0, 1), topology = (Periodic, Periodic, Bounded)); julia> k_top = Nz + 1; julia> c = ZFaceField(grid) 1×1×8 Field{Center, Center, Face} on RectilinearGrid on CPU ├── grid: 1×1×7 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×3 halo ├── boundary conditions: FieldBoundaryConditions │ └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: Nothing, top: Nothing, immersed: ZeroFlux └── data: 3×3×14 OffsetArray(::Array{Float64, 3}, 0:2, 0:2, -2:11) with eltype Float64 with indices 0:2×0:2×-2:11 └── max=0.0, min=0.0, mean=0.0 julia> set!(c, (x, y, z) -> rand()) 1×1×8 Field{Center, Center, Face} on RectilinearGrid on CPU ├── grid: 1×1×7 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×3 halo ├── boundary conditions: FieldBoundaryConditions │ └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: Nothing, top: Nothing, immersed: ZeroFlux └── data: 3×3×14 OffsetArray(::Array{Float64, 3}, 0:2, 0:2, -2:11) with eltype Float64 with indices 0:2×0:2×-2:11 └── max=0.72276, min=0.0939695, mean=0.428961 julia> cv = view(c, :, :, 1+1:k_top-1) 1×1×6 Field{Center, Center, Face} on RectilinearGrid on CPU ├── grid: 1×1×7 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×3 halo ├── boundary conditions: FieldBoundaryConditions │ └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: Nothing, top: Nothing, immersed: ZeroFlux ├── indices: (:, :, 2:7) └── data: 3×3×6 OffsetArray(view(::Array{Float64, 3}, :, :, 5:10), 0:2, 0:2, 2:7) with eltype Float64 with indices 0:2×0:2×2:7 └── max=0.72276, min=0.0939695, mean=0.503409 julia> interior(cv) 1×1×6 view(::Array{Float64, 3}, 2:2, 2:2, 5:10) with eltype Float64: [:, :, 1] = 0.6491566506900277 [:, :, 2] = 0.5118574586205709 [:, :, 3] = 0.6214555441666034 [:, :, 4] = 0.09396949649555086 [:, :, 5] = 0.7227599953730516 [:, :, 6] = 0.4212556767008805 julia> interior(cv, :, :, k_top-1) ERROR: BoundsError: attempt to access 1×1×6 view(::Array{Float64, 3}, 2:2, 2:2, 5:10) with eltype Float64 at index [1:1, 1:1, 7] Stacktrace: [1] throw_boundserror(A::SubArray{Float64, 3, Array{Float64, 3}, Tuple{UnitRange{Int64}, UnitRange{Int64}, UnitRange{Int64}}, false}, I::Tuple{Base.Slice{Base.OneTo{Int64}}, Base.Slice{Base.OneTo{Int64}}, Int64}) @ Base ./abstractarray.jl:737 [2] checkbounds @ ./abstractarray.jl:702 [inlined] [3] view(::SubArray{Float64, 3, Array{Float64, 3}, Tuple{UnitRange{Int64}, UnitRange{Int64}, UnitRange{Int64}}, false}, ::Function, ::Function, ::Int64) @ Base ./subarray.jl:184 [4] interior(::Field{Center, Center, Face, Nothing, RectilinearGrid{…}, Tuple{…}, OffsetArrays.OffsetArray{…}, Float64, FieldBoundaryConditions{…}, Nothing, Oceananigans.Fields.FieldBoundaryBuffers{…}}, ::Function, ::Function, ::Int64) @ Oceananigans.Fields /Users/Sid/Library/CloudStorage/Dropbox/StudyFolder/PostDocMITDesktop/Codes/Oceananigans/extend-parent-indices/src/Fields/field.jl:395 [5] top-level scope @ REPL[20]:1 Some type information was truncated. Use `show(err)` to see complete types.

Currently, the interior indices start from 1 and extend up to the size of the interior of the field.

I guess I'm not totally sure what the test is contributing

navidcy · 2024-05-15T08:09:25Z

test/test_field.jl

+            grid = RectilinearGrid(arch, FT, size=(Nx, Ny, Nz), x=(0, 1), y=(0, 1), z=(0, 1), topology = (Periodic, Periodic, ZTopology))
+            Hx, Hy, Hz = halo_size(grid)
+
+            k_top = (FieldType == ZFaceField && ZTopology == Bounded) ? Nz+1 : Nz 


there is a function for this; you need to do using Oceananigans.Grids: total_length and then do

Suggested change

k_top = (FieldType == ZFaceField && ZTopology == Bounded) ? Nz+1 : Nz

k_top = total_length(location(c, 3)(), topology(c, 3)(), size(grid, 3))

but after you create the field

navidcy · 2024-05-15T08:13:57Z

test/test_field.jl

+
+            # First test that the regular view is correct
+            cv = view(c, :, :, 1+1:k_top-1)
+            @test size(cv) == (1, 1, k_top-2)


Suggested change

@test size(cv) == (1, 1, k_top-2)

@test size(cv) == (Nx, Ny, k_top-2)

navidcy · 2024-05-15T08:14:12Z

test/test_field.jl

+            # First test that the regular view is correct
+            cv = view(c, :, :, 1+1:k_top-1)
+            @test size(cv) == (1, 1, k_top-2)
+            @test size(parent(cv)) == (1+2Hx, 1+2Hy, k_top-2)


Suggested change

@test size(parent(cv)) == (1+2Hx, 1+2Hy, k_top-2)

@test size(parent(cv)) == (Nx+2Hx, Ny+2Hy, k_top-2)

navidcy · 2024-05-15T08:14:35Z

test/test_field.jl

+            cv = view(c, :, :, 1+1:k_top-1)
+            @test size(cv) == (1, 1, k_top-2)
+            @test size(parent(cv)) == (1+2Hx, 1+2Hy, k_top-2)
+            CUDA.@allowscalar @test all(cv[1, 1, i] == c[1, 1, i] for i in 1+1:k_top-1)


Suggested change

CUDA.@allowscalar @test all(cv[1, 1, i] == c[1, 1, i] for i in 1+1:k_top-1)

CUDA.@allowscalar @test all(cv[i, j, k] == c[i, j, k] for k in 1+1:k_top-1, j in 1:Ny, i in 1:Nx)

@siddharthabishnu could you remove all the rest hardcoded 1's below in the tests? (unless they are not hard coded and I missed it...)

navidcy

made few suggestions

siddharthabishnu · 2024-05-15T23:54:43Z

made few suggestions

Incorporated in commit d71e14a.

glwagner · 2024-05-16T04:15:31Z

good to merge by me

navidcy · 2024-05-16T07:19:57Z

I'm merging

Just a hack for now!

3a0cc14

This was referenced May 1, 2024

Bug in determination of indices of parent array by parent_index_range #3572

Closed

Employ new halo-filling functions to update halos of prognostic variables on the cubed sphere #3570

Draft

Momentum dynamics validation on the CubedSphere #3306

Open

siddharthabishnu changed the title ~~Extend parent array to variables without z-halos~~ Extend Parent Array to Support 1D or 2D Fields Without Halos in Nonexistent Dimensions May 1, 2024

glwagner reviewed May 1, 2024

View reviewed changes

src/Fields/field.jl Outdated Show resolved Hide resolved

siddharthabishnu and others added 2 commits May 1, 2024 13:15

A real solution this time!

bbfc656

Merge branch 'main' into sb/extend-parent-indices

4d66f73

glwagner changed the title ~~Extend Parent Array to Support 1D or 2D Fields Without Halos in Nonexistent Dimensions~~ Bugfix for parent_index_range May 1, 2024

glwagner marked this pull request as ready for review May 1, 2024 17:29

navidcy added the bug 🐞 Even a perfect program still has bugs label May 2, 2024

Make views of views work

f217b80

Merge branch 'main' into sb/extend-parent-indices

605a613

Create tests for views of views

e2be08b

glwagner reviewed May 4, 2024

View reviewed changes

test/test_field.jl Show resolved Hide resolved

navidcy and others added 3 commits May 8, 2024 09:27

Merge branch 'main' into sb/extend-parent-indices

20e3bff

Refine code

acc9e60

Use CUDA.@allowscalar

9543e05

siddharthabishnu added 2 commits May 9, 2024 16:55

Slightly compact code

cbc18d0

Merge 'origin/main' into extend-parent-indices

58cd98e

glwagner reviewed May 13, 2024

View reviewed changes

siddharthabishnu and others added 3 commits May 14, 2024 22:02

Delete redundant comments

7d68c02

Use CUDA.@allowscalar only where necessary

427cb51

Merge branch 'main' into sb/extend-parent-indices

6dc5363

navidcy reviewed May 15, 2024

View reviewed changes

Remove hardcoded 1s and refine code

d71e14a

glwagner approved these changes May 16, 2024

View reviewed changes

navidcy merged commit 11c317d into main May 16, 2024
46 checks passed

navidcy deleted the sb/extend-parent-indices branch May 16, 2024 07:20

	@test_throws BoundsError cvvv[:, :, k_top-2:k_top]
	@test_throws BoundsError cvvv[:, :, k_top-2:k_top]

	k_top = (FieldType == ZFaceField && ZTopology == Bounded) ? Nz+1 : Nz
	k_top = total_length(location(c, 3)(), topology(c, 3)(), size(grid, 3))

	@test size(cv) == (1, 1, k_top-2)
	@test size(cv) == (Nx, Ny, k_top-2)

	@test size(parent(cv)) == (1+2Hx, 1+2Hy, k_top-2)
	@test size(parent(cv)) == (Nx+2Hx, Ny+2Hy, k_top-2)

	CUDA.@allowscalar @test all(cv[1, 1, i] == c[1, 1, i] for i in 1+1:k_top-1)
	CUDA.@allowscalar @test all(cv[i, j, k] == c[i, j, k] for k in 1+1:k_top-1, j in 1:Ny, i in 1:Nx)

Bugfix for parent_index_range #3573

Bugfix for parent_index_range #3573

Conversation

siddharthabishnu commented May 1, 2024

glwagner commented May 1, 2024

glwagner commented May 1, 2024

navidcy commented May 2, 2024

glwagner commented May 2, 2024

glwagner commented May 2, 2024 • edited

glwagner commented May 2, 2024

glwagner commented May 2, 2024 • edited

navidcy commented May 2, 2024

glwagner commented May 2, 2024

glwagner commented May 2, 2024

glwagner commented May 2, 2024

glwagner commented May 2, 2024

navidcy commented May 2, 2024

navidcy commented May 2, 2024

navidcy commented May 2, 2024 • edited

glwagner commented May 2, 2024

navidcy commented May 2, 2024

glwagner May 4, 2024

Choose a reason for hiding this comment

glwagner commented May 4, 2024

navidcy commented May 8, 2024

siddharthabishnu commented May 9, 2024 • edited

siddharthabishnu commented May 9, 2024

siddharthabishnu commented May 13, 2024

glwagner May 13, 2024

Choose a reason for hiding this comment

siddharthabishnu May 15, 2024 • edited

Choose a reason for hiding this comment

glwagner May 13, 2024

Choose a reason for hiding this comment

siddharthabishnu May 15, 2024 • edited

Choose a reason for hiding this comment

glwagner May 16, 2024

Choose a reason for hiding this comment

navidcy May 15, 2024 • edited

Choose a reason for hiding this comment

navidcy May 15, 2024

Choose a reason for hiding this comment

navidcy May 15, 2024

Choose a reason for hiding this comment

navidcy May 15, 2024 • edited

Choose a reason for hiding this comment

navidcy May 15, 2024 • edited

Choose a reason for hiding this comment

navidcy left a comment

Choose a reason for hiding this comment

siddharthabishnu commented May 15, 2024 • edited

glwagner commented May 16, 2024

navidcy commented May 16, 2024

Bugfix for `parent_index_range` #3573

Bugfix for `parent_index_range` #3573

glwagner commented May 2, 2024 •

edited

glwagner commented May 2, 2024 •

edited

navidcy commented May 2, 2024 •

edited

siddharthabishnu commented May 9, 2024 •

edited

siddharthabishnu May 15, 2024 •

edited

siddharthabishnu May 15, 2024 •

edited

navidcy May 15, 2024 •

edited

navidcy May 15, 2024 •

edited

navidcy May 15, 2024 •

edited

siddharthabishnu commented May 15, 2024 •

edited