NamedTuple-Field r-operations are broken

[charleskawczynski]

MWE:

import ClimaCore.Geometry as Geometry
import ClimaCore.Domains as Domains
import ClimaCore.Meshes as Meshes
import ClimaCore.Topologies as Topologies
import ClimaCore.Spaces as Spaces
import ClimaCore.Fields as Fields
import ClimaCore.Operators as Operators

function foo(::Type{FT}) where {FT}
    nelems = 20
    npoly = 3
    domain = Domains.SphereDomain(FT(6371220))
    mesh = Meshes.EquiangularCubedSphere(domain, nelems)
    topology = Topologies.Topology2D(mesh)
    quad = Spaces.Quadratures.GLL{npoly + 1}()
    space = Spaces.SpectralElementSpace2D(topology, quad)
    coords = Fields.coordinate_field(space)

    ρ = ones(space)
    q = map(coords) do coord
        (; gaussian_hills=FT(0), cosine_bells = FT(0), slotted_cylinders = FT(0))
    end
    init_state = Fields.FieldVector(; ρ, ρq = ρ .* q)

    # current wind vector (Section 2.3)
    current_wind_vector = Fields.Field(Geometry.UVVector{FT}, space)
    wind_vector(coord, t) = Geometry.UVVector(FT(0), FT(0))

    div = Operators.Divergence()
    function T_lim!(tendency, state, _, t)
        @. current_wind_vector = wind_vector(coords, t)
        @. tendency.ρ = -div(state.ρ * current_wind_vector)
        @. tendency.ρq = -div(state.ρq * current_wind_vector)
        return nothing
    end
    tendency = similar(init_state)
    state = init_state
    T_lim!(tendency, state, (;), FT(0))
end
foo(Float64)
nothing

This used to work in version 0.10.39 (found in ClimaTimesteppers).

[charleskawczynski]

@simonbyrne, it may be that this limiter test in ClimaTimeSteppers was silently not failing (until recently). Can you confirm?

In particular:

@. tendency.ρq = -div(state.ρq * current_wind_vector)

looks to me like an ill-defined operation (state.ρq is a NamedTuple-valued Field). However, I suppose that this could just operate point-wise over the values in the NamedTuple.

[charleskawczynski]

Stack trace is:

julia> include("../cc_inference_problem.jl")
ERROR: LoadError: MethodError: rmaptype(::typeof(ClimaCore.Geometry.divergence_result_type), ::Type{Union{}}) is ambiguous.

Candidates:
  rmaptype(fn::F, ::Type{T}) where {F, T<:Tuple}
    @ ClimaCore.RecursiveApply dev/ClimaCore.jl/src/RecursiveApply/RecursiveApply.jl:74
  rmaptype(fn::F, ::Type{T}) where {F, names, Tup, T<:NamedTuple{names, Tup}}
    @ ClimaCore.RecursiveApply dev/ClimaCore.jl/src/RecursiveApply/RecursiveApply.jl:76

Possible fix, define
  rmaptype(::F, ::Type{Union{}}) where F

Stacktrace:
  [1] operator_return_eltype(op::ClimaCore.Operators.Divergence{(1, 2)}, #unused#::Type{Union{}})
    @ ClimaCore.Operators dev/ClimaCore.jl/src/Operators/spectralelement.jl:709
  [2] apply_operator
    @ dev/ClimaCore.jl/src/Operators/spectralelement.jl:753 [inlined]
  [3] resolve_operator
    @ dev/ClimaCore.jl/src/Operators/spectralelement.jl:218 [inlined]
  [4] _resolve_operator_args
    @ dev/ClimaCore.jl/src/Operators/spectralelement.jl:235 [inlined]
  [5] resolve_operator
    @ dev/ClimaCore.jl/src/Operators/spectralelement.jl:224 [inlined]
  [6] copyto_slab!
    @ dev/ClimaCore.jl/src/Operators/spectralelement.jl:193 [inlined]
  [7] #1
    @ dev/ClimaCore.jl/src/Operators/spectralelement.jl:178 [inlined]
  [8] byslab
    @ dev/ClimaCore.jl/src/Fields/indices.jl:218 [inlined]
  [9] byslab
    @ dev/ClimaCore.jl/src/Fields/indices.jl:208 [inlined]
 [10] copyto!
    @ dev/ClimaCore.jl/src/Operators/spectralelement.jl:176 [inlined]
 [11] materialize!
    @ ./broadcast.jl:884 [inlined]
 [12] materialize!
    @ ./broadcast.jl:881 [inlined]
 [13] (::var"#T_lim!#65"{ClimaCore.Operators.Divergence{()}, var"#wind_vector#64"{Float64, Float64, Floa
...
Next line points to: `@. tendency.ρq = -div(state.ρq * current_wind_vector)`

[simonbyrne]

looks to me like an ill-defined operation (state.ρq is a NamedTuple-valued Field). However, I suppose that this could just operate point-wise over the values in the NamedTuple.

yes, that's what it does

[simonbyrne]

ah, maybe the inference broke for that. We don't make use of that at the moment.

[charleskawczynski]

This is actually a more general issue that our RecursiveApply operations are broken on NamedTuple-valued fields. MWE:

using Revise
import ClimaCore.Geometry as Geometry
import ClimaCore.Domains as Domains
import ClimaCore.Meshes as Meshes
import ClimaCore.Topologies as Topologies
import ClimaCore.Spaces as Spaces
import ClimaCore.Fields as Fields
import ClimaCore.Operators as Operators

function get_args(::Type{FT}) where {FT}
    nelems = 20
    npoly = 3
    domain = Domains.SphereDomain(FT(6371220))
    mesh = Meshes.EquiangularCubedSphere(domain, nelems)
    topology = Topologies.Topology2D(mesh)
    quad = Spaces.Quadratures.GLL{npoly + 1}()
    space = Spaces.SpectralElementSpace2D(topology, quad)
    coords = Fields.coordinate_field(space)

    ρ = ones(space)
    q = map(coords) do coord
        (; gaussian_hills=FT(0), cosine_bells = FT(0), slotted_cylinders = FT(0))
    end
    init_state = Fields.FieldVector(; ρ, ρq = ρ .* q)

    # current wind vector (Section 2.3)
    current_wind_vector = Fields.Field(Geometry.UVVector{FT}, space)
    p = (;current_wind_vector)
    state = init_state
    return (; state, p)
end
(;state, p) = get_args(Float64)
(;current_wind_vector) = p
import ClimaCore.RecursiveApply: ⊠
ρq = state.ρq
Y = ρq .⊠ current_wind_vector

nothing

[charleskawczynski]

A simpler MWE still:

using Test
import ClimaCore.Geometry as Geometry
import ClimaCore.RecursiveApply as RecursiveApply
FT = Float64
nt = (; a = FT(1), b = FT(2))
uv = Geometry.UVVector(FT(1), FT(2))
@test rz = RecursiveApply.rmap(*, nt, uv)
@test typeof(rz) == NamedTuple{(:a, :b), Tuple{UVVector{FT}, UVVector{FT}}}
@test @inferred RecursiveApply.rmap(*, nt, uv)
@test rz.a.u == 1
@test rz.a.v == 2
@test rz.b.u == 1
@test rz.b.v == 4

[charleskawczynski]

closed by #1454