CliMA / ClimaCore.jl

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

Clarify 2d spaces #1578

Open LenkaNovak opened 8 months ago

LenkaNovak commented 8 months ago

Is your feature request related to a problem? Please describe. Currently, there are two types of 2d fields

when queried, both print the same info, e.g.:

SpectralElementSpace2D:
  context: SingletonCommsContext using CPUSingleThreaded
  mesh: 4×4×6-element EquiangularCubedSphere of SphereDomain: radius = 6.371e6
  quadrature: 4-point Gauss-Legendre-Lobatto quadrature

but 2. still contains it's original 3d info and jacobian, etc. This may cause problems when dispatching on space type and/or calculating area integrals using sum.

Describe the solution you'd like It would be good to make 2. a 3d space and allow an additional 2d reduction from that.

Min demo example

using ClimaCore: Fields, Domains, Meshes, Topologies, Spaces, Geometry, Utilities
function create_space(
    FT;
    comms_ctx = ClimaComms.SingletonCommsContext(),
    R = FT(6371e3),
    ne = 4,
    polynomial_degree = 3,
    nz = 1,
)
    domain = Domains.SphereDomain(R)
    mesh = Meshes.EquiangularCubedSphere(domain, ne)

    if comms_ctx isa ClimaComms.SingletonCommsContext
        topology = Topologies.Topology2D(mesh, Topologies.spacefillingcurve(mesh))
    else
        topology = Topologies.DistributedTopology2D(comms_ctx, mesh, Topologies.spacefillingcurve(mesh))
    end

    Nq = polynomial_degree + 1
    quad = Spaces.Quadratures.GLL{Nq}()
    sphere_space = Spaces.SpectralElementSpace2D(topology, quad)

    if nz > 1
        vertdomain =
            Domains.IntervalDomain(Geometry.ZPoint{FT}(0), Geometry.ZPoint{FT}(100); boundary_tags = (:bottom, :top))
        vertmesh = Meshes.IntervalMesh(vertdomain, nelems = nz)
        vert_center_space = Spaces.CenterFiniteDifferenceSpace(vertmesh)
        return Spaces.ExtrudedFiniteDifferenceSpace(sphere_space, vert_center_space)
    else
        return sphere_space
    end
end

FT = Float64
space_2d = create_space(FT, nz = 1, R = 6371e3)
space_3d = create_space(FT, nz = 3, R = 6371e3)
space_2d_from_3d = Fields.level(space_3d, 1)

space_2d isa Spaces.SpectralElementSpace2D # true
space_2d_from_3d isa Spaces.SpectralElementSpace2D # true

expected_sum = 4π * 6371e3^2 # 5.1006447190978825e14
sum(ones(space_2d)) # 5.100617440612155e14
sum(ones(space_2d_from_3d)) # 1.700205813537382e16
LenkaNovak commented 8 months ago

Currently, to get the expected surface integrals, the following corrections needs to be applied, depending on whether the original field was center or face:

sum(ones(space_2d_from_3d) ./ Fields.dz_field(space_2d_from_3d) ) # 5.100617440612155e14 (if derived from CenterFiniteDifferenceSpace)

sum(ones(space_2d_from_3d) ./ Fields.dz_field(space_2d_from_3d) .* 2 ) # 5.100617440612155e14 (if derived from FaceFiniteDifferenceSpace)

Whilst we can write convenience functions for this, it seems that there could be a more sophisticated approach that distinguishes between the 2d/3d spaces more clearly.