Unintuitive behavior of sum

[kmdeck]

Describe the bug

I am trying to use sum to compute integrals of different fields over different spaces. In particular, I am taking sum(ones(space)) - which should be volume or area or column or "point" integrals. But,

• I get funny results when using level of the 3d space vs the horizontal space for the spherical shell. For example, I would expect that \int (ones(horizontal_space)) = 4pi R^2, and that \int ones(level) = 4pi(R+z)^2, where z is the height of the level. What \int ones(level) returns - for any level, at any height - is 4\pi R^2 dz/2, where dz is the discretization of the vertical.
• I would expect sum(point_field) = 0, but I get a nonzero result.
• I also found that using by column to take a column integral over each column of the spherical shell - returning something that lives on the horizontal space - also seemed to give unintuitive results.

To Reproduce

MWE below

using ClimaCore
using Test
@testset "Spherical Shell" begin
    FT = Float32
    radius = FT(10.0)
    height =FT(2.0)
    nelements = (10,10)
    npolynomial = 3
    vertdomain = ClimaCore.Domains.IntervalDomain(
        ClimaCore.Geometry.ZPoint(FT(0)),
        ClimaCore.Geometry.ZPoint(FT(height));
        boundary_tags = (:bottom, :top),
    )
    
    vertmesh = ClimaCore.Meshes.IntervalMesh(
        vertdomain,
        ClimaCore.Meshes.Uniform(),
        nelems = nelements[2],
    )
    vert_center_space = ClimaCore.Spaces.CenterFiniteDifferenceSpace(vertmesh)

    horzdomain = ClimaCore.Domains.SphereDomain(radius)
    horzmesh = ClimaCore.Meshes.EquiangularCubedSphere(horzdomain, nelements[1])
    horztopology = ClimaCore.Topologies.Topology2D(horzmesh)
    quad = ClimaCore.Spaces.Quadratures.GLL{npolynomial + 1}()
    horzspace = ClimaCore.Spaces.SpectralElementSpace2D(horztopology, quad)

    hv_center_space = ClimaCore.Spaces.ExtrudedFiniteDifferenceSpace(
        horzspace,
        vert_center_space,
    )
    hv_face_space = ClimaCore.Spaces.FaceExtrudedFiniteDifferenceSpace(hv_center_space)
    # Sum behavior
    ones_volume = ones(hv_center_space)
    ones_horizontal_space = ones(hv_center_space.horizontal_space)
    ones_bottom_space_level = ones(ClimaCore.Spaces.level(hv_face_space, ClimaCore.Utilities.PlusHalf(0)))
    ones_top_space_level = ones(ClimaCore.Spaces.level(hv_face_space, ClimaCore.Utilities.PlusHalf(10)))
    # Passing tests
    @test sum(ones_volume) ≈ FT(4.0*π)*radius^2.0*height 
    @test sum(ones_horizontal_space) ≈  FT(4.0*π)*radius^2.0

    # Not passing
    @test sum(ones_bottom_space_level) ≈ FT(4.0*π)*radius^2.0 
    @test sum(ones_top_space_level) ≈ FT(4.0*π)*(radius+height)^2.0
    column_integrals = zeros(horzspace)
    ClimaCore.Fields.bycolumn(horzspace) do colidx
        column_integrals[colidx] .= sum(ones_volume[colidx])
    end
    # Expect this should be ∫dz for each element? = height?
    @test sum(parent(column_integrals .- height))≈ FT(0.0)
end

@testset "Column" begin
    zlim = Float32.((0.0,10.0))
    nelements = 20
    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)
    face_space = ClimaCore.Spaces.FaceFiniteDifferenceSpace(mesh)
    bottom_space = ClimaCore.Spaces.level(face_space, ClimaCore.Utilities.PlusHalf(0))
    
    ones_column = ones(center_space)
    ones_point_bottom = ones(bottom_space)
    #Passes
    @test sum(ones_column) ≈ FT(10)
    # Does not pass
    @test sum(ones_point_bottom) ≈ FT(0) #∫dx over a point is zero
end

System details

Any relevant system information:

• Julia version 1.7
• operating system Mac

[simonbyrne]

This one has bothered me as well: I admit I don't have a good solution.

[kmdeck]

Is sum actually doing an integral over the domain? Maybe it would be clearer if there was a formula showing what sum was computing, and then a user could see how it reduces in some case to the integral and some cases not? Im not sure. I dont know what it is doing under the hood.

[valeriabarra]

@kmdeck yes, my understanding is that sum has to use the Jacobian information of the underlying domain of the field you are trying to compute the integral of.
In essence, think about a change in coordinate transformation. Say, for instance, that you switch from Cartesian to polar coordinates, then any integral you compute has to account for the change of coordinates (and that is, in fact, the Jacobian that comes out of the chain rule). This is why we have to account for WJ and we tried to make this clear in the docstring (the W is due to the quadrature point weights needed to compute integrals, since we are using spectral spaces. These are just equal to 1 for Finite Difference spaces).

So to address your second bullet point, yes, we should expect sum of a point-field to return just zero.

For the first bullet point, is it due to the fact that we use the "shallow-atmosphere" approximation for our 3D sphere? See for instance this comment in this test.

[kmdeck]

Hi @valeriabarra - thank you! It does make sense that we need the Jacobian in the integral. That said, the results from sum still dont make sense in the spherical case (to me!).

I dont think that the issue in the first test/bullet point is because of the shallow-atmosphere approx. The sum I was showing in the test is over a 2D spectral element space, so the height is only involved in that it affects the radius of the sphere we are taking an area integral of. In the shallow case, we would expect 4pi(R+h)^2 \approx 4piR^2(1+2h/R), also, and not 4piR^2 dz/2...

as a comment - making the shallow approximation makes a lot of sense for the globe, but there is no restriction on what the user can pass in for h and r, so that assumption shouldnt be used in the formula for sum anywhere...

(Also note the funny behavior with bycolumn - seperate issue)

[charleskawczynski]

FWIW, I would expect the integral over a point-space to be non-zero (in general), but it should be decreasing as the resolution increases.

[charleskawczynski]

One reason why I think it's somewhat dangerous/unintuitive to define sum as an integral is that a user-friendly API of computing integrals (or any weighted sum) only make sense on fields, and not data-layouts (where users are responsible for properly incorporating metric terms).

[charleskawczynski]

Here is an updated reproducer:

using ClimaCore.CommonSpaces
using ClimaCore: Spaces
using ClimaCore: Spaces
using ClimaCore: Utilities
using ClimaCore: Fields
using ClimaCore
using Test
@testset "Spherical Shell" begin
    FT = Float32
    radius = FT(10.0)
    height = FT(2.0)
    n_quad_points = 3

    hv_center_space = ExtrudedCubedSphereSpace(FT;
        z_elem = 10,
        z_min = 0,
        z_max = height,
        radius = 10,
        h_elem = 10,
        n_quad_points = 4,
        staggering = CellCenter()
    )
    hv_face_space = Spaces.face_space(hv_center_space)
    # Sum behavior
    ones_volume = ones(hv_center_space)
    hspace = Spaces.horizontal_space(hv_center_space)
    ones_horizontal_space = ones(Spaces.horizontal_space(hv_center_space))
    ones_bottom_space_level = ones(Spaces.level(hv_face_space, Utilities.PlusHalf(0)))
    ones_top_space_level = ones(Spaces.level(hv_face_space, Utilities.PlusHalf(10)))
    # Passing tests
    @test sum(ones_volume) ≈ FT(4.0*π)*radius^2.0*height
    @test sum(ones_horizontal_space) ≈  FT(4.0*π)*radius^2.0

    # Not passing
    @test sum(ones_bottom_space_level) ≈ FT(4.0*π)*radius^2.0
    @test sum(ones_top_space_level) ≈ FT(4.0*π)*(radius+height)^2.0
    column_integrals = zeros(hspace)
    Fields.bycolumn(hspace) do colidx
        column_integrals[colidx] .= sum(ones_volume[colidx])
    end
    # Expect this should be ∫dz for each element? = height?
    @test sum(parent(column_integrals .- height))≈ FT(0.0)
end

@testset "Column" begin
    FT = Float32
    center_space = ColumnSpace(FT;
        z_elem = 20,
        z_min = 0,
        z_max = 10,
        staggering = CellCenter()
    )
    face_space = Spaces.face_space(center_space)
    bottom_space = Spaces.level(face_space, Utilities.PlusHalf(0))

    ones_column = ones(center_space)
    ones_point_bottom = ones(bottom_space)
    #Passes
    @test sum(ones_column) ≈ FT(10)
    # Does not pass
    @test sum(ones_point_bottom) ≈ FT(0) #∫dx over a point is zero
end