CumulativeIntegral throws error when integrated over a reduced Field

[tomchor]

It appears there's a bug with CumulativeIntegral when used on reduced Fields.

MWE:

julia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));

julia> c = CenterField(grid);

julia> c̄ = Field(Average(c, dims=(1, 2)));

julia> ∫c̄ = CumulativeIntegral(c̄, dims=3)
ERROR: UndefVarError: `Δzᵃᵃᶜ` not defined
Stacktrace:
 [1] metric_function
   @ /glade/work/tomasc/.julia/packages/Oceananigans/Dj1rO/src/AbstractOperations/grid_metrics.jl:31 [inlined]
 [2] Oceananigans.AbstractOperations.GridMetricOperation(L::Tuple{…}, metric::Function, grid::RectilinearGrid{…})
   @ Oceananigans.AbstractOperations /glade/work/tomasc/.julia/packages/Oceananigans/Dj1rO/src/AbstractOperations/grid_metrics.jl:78
 [3] *(Lc::Tuple{…}, a::Field{…}, m::typeof(Oceananigans.Operators.Δz))
   @ Oceananigans.AbstractOperations /glade/work/tomasc/.julia/packages/Oceananigans/Dj1rO/src/AbstractOperations/binary_operations.jl:122
 [4] *
   @ /glade/work/tomasc/.julia/packages/Oceananigans/Dj1rO/src/AbstractOperations/binary_operations.jl:126 [inlined]
 [5] (Accumulation{…})(field::Field{…}; dims::Int64, reverse::Bool, condition::Nothing, mask::Int64)
   @ Oceananigans.AbstractOperations /glade/work/tomasc/.julia/packages/Oceananigans/Dj1rO/src/AbstractOperations/metric_field_reductions.jl:170
 [6] top-level scope
   @ REPL[20]:1
Some type information was truncated. Use `show(err)` to see complete types.

The same thing happens when using Integral instead of Average.

[glwagner]

A little more direct MWE:

julia> c = Field{Nothing, Nothing, Center}(grid);

julia> ∫c = CumulativeIntegral(c, dims=3)
ERROR: UndefVarError: `Δzᵃᵃᶜ` not defined
Stacktrace:
 [1] metric_function
   @ ~/Projects/dev/Oceananigans.jl/src/AbstractOperations/grid_metrics.jl:31 [inlined]
 [2] Oceananigans.AbstractOperations.GridMetricOperation(L::Tuple{…}, metric::Function, grid::RectilinearGrid{…})
   @ Oceananigans.AbstractOperations ~/Projects/dev/Oceananigans.jl/src/AbstractOperations/grid_metrics.jl:78
 [3] *(Lc::Tuple{…}, a::Field{…}, m::typeof(Oceananigans.Operators.Δz))
   @ Oceananigans.AbstractOperations ~/Projects/dev/Oceananigans.jl/src/AbstractOperations/binary_operations.jl:122
 [4] *
   @ ~/Projects/dev/Oceananigans.jl/src/AbstractOperations/binary_operations.jl:126 [inlined]
 [5] (Accumulation{…})(field::Field{…}; dims::Int64, reverse::Bool, condition::Nothing, mask::Int64)
   @ Oceananigans.AbstractOperations ~/Projects/dev/Oceananigans.jl/src/AbstractOperations/metric_field_reductions.jl:170
 [6] top-level scope
   @ REPL[73]:1
Some type information was truncated. Use `show(err)` to see complete types.

The issue I think is that because c has no location in x, y, we are looking for a metric that also has no location in x, y for the vertical integral. However, no such metric exists in general (for example if we are on an ImmersedBoundaryGrid with PartialCellBottom). What is the behavior that we would like to implement, in the case that Δz depends on x, y? We have to decide on an effective vertical metric for a reduced field.

[tomchor]

Yes, although for non immersed grids we don't have that problem and this should work, no?

Also, a quick note. While CumulativeIntegral(Field(Average(c, dims=(1, 2))), dims=3) fails, using CumulativeIntegral first and then averaging in the horizontal works:

julia> Average(Field(CumulativeIntegral(c, dims=3)), dims=(1, 2))
Average of 4×4×4 Field{Center, Center, Center} on RectilinearGrid on CPU over dims (1, 2)
└── operand: 4×4×4 Field{Center, Center, Center} on RectilinearGrid on CPU
    └── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo

[glwagner]

Yes I think if you want to make some special cases for grids work then you can extend this constructor?

Oceananigans.jl/src/AbstractOperations/grid_metrics.jl

Line 78 in a8cf967

GridMetricOperation(L, metric, grid) = GridMetricOperation{L[1], L[2], L[3]}(metric_function(L, metric), grid)

[glwagner]

Yes, although for non immersed grids we don't have that problem and this should work, no?

Also, a quick note. While CumulativeIntegral(Field(Average(c, dims=(1, 2))), dims=3) fails, using CumulativeIntegral first and then averaging in the horizontal works:

julia> Average(Field(CumulativeIntegral(c, dims=3)), dims=(1, 2))
Average of 4×4×4 Field{Center, Center, Center} on RectilinearGrid on CPU over dims (1, 2)
└── operand: 4×4×4 Field{Center, Center, Center} on RectilinearGrid on CPU
└── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo

Is this correct? I was playing around with CumulativeIntegral and it seems wrong that the location of the integral is the same as the location of the Field. I expected the location to flip, eg the same way it flips for a derivative. With that logic the cumulative integral of c would be at Face in the vertical.