CliMA / ClimaCore.jl

CliMA model dycore
https://clima.github.io/ClimaCore.jl/dev
Apache License 2.0
87 stars 8 forks source link

Strange behavior of `similar` and broadcast with Fields of IJFH datalayout #598

Open jb-mackay opened 2 years ago

jb-mackay commented 2 years ago

Was playing around with Field broadcasting and stumbled into some unexpected behavior. My setup (from test/Operators/remapping.jl) is below.

field1 = Float64-valued Field:
  [1.0, 1.0, 1.0, 1.0]

field2 = Float64-valued Field:
  [2.0, 2.0, 2.0, 2.0]

When I write field2 .= field1 I expect to have

field2 = Float64-valued Field:
  [1.0, 1.0, 1.0, 1.0]

but instead I have

field2 = Float64-valued Field:
  [1.0, 2.0, 2.0, 2.0]

where only the first nodal value is changed. Furthermore, similar(field1) returns a Float64-valued Field: [2.17214e-314].

So again, only a single nodal value is represented.

To reproduce:

using ClimaCore:
    Domains, Meshes, Topologies, Geometry, Operators, Spaces, Fields
using ClimaCore.Operators: local_weights, LinearRemap, remap, remap!
using ClimaCore.Spaces: Quadratures
using ClimaCore.DataLayouts: IJFH
using IntervalSets

FT = Float64

function make_space(domain::Domains.IntervalDomain, nq, nelems = 1)
    nq == 1 ? (quad = Quadratures.GL{1}()) : (quad = Quadratures.GLL{nq}())
    mesh = Meshes.IntervalMesh(domain; nelems = nelems)
    topo = Topologies.IntervalTopology(mesh)
    space = Spaces.SpectralElementSpace1D(topo, quad)
    return space
end

domain = Domains.IntervalDomain(
                Geometry.XPoint(0.0) .. Geometry.XPoint(1.0),
                boundary_tags = (:left, :right),
            )

source = make_space(domain, 1, 4)
n = length(source.local_geometry)
field1 = Fields.Field(IJFH{FT, 1}(ones(1, 1, n, 1)), source)
field2 = Fields.Field(IJFH{FT, 1}(2 * ones(1, 1, n, 1)), source)

field2 .= field1
field3 = similar(field1)

I'll note that I should update to using Fields.ones(source) going forward, and using that this issue does not come up. So there is something with how I construct these fields that is causing this behavior.

jakebolewski commented 2 years ago

similar is for AbstractArray subtypes and Field's are explicitly not arrays, so I suspect the generic version treats the Field as a scalar by default but we overload enough methods of AbstractArray interface that you get some result (the backing data layouts have length 1). Assigning fields to each other on the same space is a bug but you would have to special case that and not use the generic broadcast fallback.

jb-mackay commented 2 years ago

That's a good point about the data layout length.

Assigning fields to each other on the same space is a bug but you would have to special case that and not use the generic broadcast fallback.

What methods need special casing here?

In general, what are the to-dos that arise from this?