CliMA / CliMA/ClimaCore.jl

Unintuitive behavior of `sum`

Open
#943 7 comments 0 reactions 0 assignees View on GitHub
bug
Dominant language
Julia
Stars
117
Forks
19
Avg merge
3d 4h
Merged PRs (30d)
41

Description

## 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
```julia
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

Contributor guide

No contributing guide indexed for this repository

Assessment

This issue has not been assessed yet.

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.