ClimaCore.jl icon indicating copy to clipboard operation
ClimaCore.jl copied to clipboard

Document how to interpolate when indexing into tuple fields

Open kmdeck opened this issue 1 year ago • 2 comments

Is your feature request related to a problem? Please describe. It's unclear what the rules are for indexing into a Field of Tuples. For example,

using ClimaCore
zlim = (1,2)
nelements = 10
FT = Float32
boundary_tags = (:bottom, :top)
column = ClimaCore.Domains.IntervalDomain(
    ClimaCore.Geometry.ZPoint{FT}(zlim[1]),
    ClimaCore.Geometry.ZPoint{FT}(zlim[2]);
    boundary_tags = boundary_tags,
)
mesh = ClimaCore.Meshes.IntervalMesh(column; nelems = nelements)
center_space = ClimaCore.Spaces.CenterFiniteDifferenceSpace(mesh)
coords = ClimaCore.Fields.coordinate_field(center_space)
types = (NTuple{10, FT},)
fields = map(types) do (T)
    zero_instance = ClimaCore.RecursiveApply.rzero(T)
    map(_ -> zero_instance, coords)
end
field = fields[1]
f(x) = x+x

## Indexing outside of loops

field.:1 .= FT(1) # runs
id = 2+3
field.:($id) .= FT(3) # runs
field.:($(2+3)) # does not run
field.:($(2+3)) .= FT(0) # runs
@. field.:($(2+3)) = FT(0) # does not run
field.:($$(2+3)) # does not run
@. field.:($$(2+3)) = FT(0) # runs

@. field.:1 = f(field.:1) # runs
@. field.:($$(2+3)) = f(field.:1) # runs
@. field.:($$(2+3)) = f(field.:($$(2+3))) # does not run

The rules are slightly different in loops: This all runs

for i in 1:9
    ip1 = i+1
    field.:($i) .= field.:($ip1)
    @. field.:($$i) = FT(0)
    @. field.:($$ip1) = FT(0)
    @. field.:($$ip1) =  field.:($$ip1)
    @. field.:($$(i+1)) = field.:($$ip1)
end

This does not

for i in 1:9
    @. field.:($$(i+1)) = field.:($$(i+1))
end

Describe the solution you'd like Documentation and explanation in a central repository, since this type of indexing is done in multiple repos.

kmdeck avatar Sep 07 '23 17:09 kmdeck

@juliasloan25 @charleskawczynski

kmdeck avatar Sep 07 '23 17:09 kmdeck

One general rule that I think I'm seeing is that if we interpolate integers into a getproperty call (.:($i)) is that we need an extra interpolation (extra $) if it's inside a broadcast expression.

charleskawczynski avatar Oct 07 '23 02:10 charleskawczynski