CliMA / ClimaCore.jl

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

Not all operators take fields as boundary conditions (low priority) #886

Open kmdeck opened 2 years ago

kmdeck commented 2 years ago

Is your feature request related to a problem? Please describe. The divf2c operator (at least) takes a field (on the boundary space) as a boundary condition:

using ClimaCore

FT = Float64
zmin = FT(1.0)
zmax = FT(2.0)
xlim = FT.((0.0, 10.0))
ylim = FT.((0.0, 1.0))
zlim = FT.((zmin, zmax))
nelements = (1, 1, 5)
npolynomial = 3
domain_x = ClimaCore.Domains.IntervalDomain(
    ClimaCore.Geometry.XPoint(xlim[1]),
    ClimaCore.Geometry.XPoint(xlim[2]);
    periodic = true,
)
domain_y = ClimaCore.Domains.IntervalDomain(
    ClimaCore.Geometry.YPoint(ylim[1]),
    ClimaCore.Geometry.YPoint(ylim[2]);
    periodic = true,
)
plane = ClimaCore.Domains.RectangleDomain(domain_x, domain_y)

mesh =
    ClimaCore.Meshes.RectilinearMesh(plane, nelements[1], nelements[2])
grid_topology = ClimaCore.Topologies.Topology2D(mesh)
quad = ClimaCore.Spaces.Quadratures.GLL{npolynomial + 1}()
horzspace = ClimaCore.Spaces.SpectralElementSpace2D(grid_topology, quad)

vertdomain = ClimaCore.Domains.IntervalDomain(
    ClimaCore.Geometry.ZPoint(zlim[1]),
    ClimaCore.Geometry.ZPoint(zlim[2]);
    boundary_tags = (:bottom, :top),
)
vertmesh = ClimaCore.Meshes.IntervalMesh(vertdomain, nelems = nelements[3])
vert_center_space = ClimaCore.Spaces.CenterFiniteDifferenceSpace(vertmesh)

hv_center_space =
    ClimaCore.Spaces.ExtrudedFiniteDifferenceSpace(horzspace, vert_center_space)

surface_field = ClimaCore.Fields.zeros(horzspace)

ψ = ClimaCore.Fields.ones(hv_center_space)
value = surface_field # same type/instance of underlying space as horizontal space of \psi
gradc2f_no_bc = ClimaCore.Operators.GradientC2F()
divf2c = ClimaCore.Operators.DivergenceF2C(
                  top = ClimaCore.Operators.SetValue(ClimaCore.Geometry.WVector.(value)),
                  bottom = ClimaCore.Operators.SetValue(ClimaCore.Geometry.WVector(FT(0.0))),
              )
@. divf2c(gradc2f_no_bc(ψ)) # runs

This does not work for gradc2f:

gradc2f = ClimaCore.Operators.GradientC2F(top = ClimaCore.Operators.SetValue(value), bottom = ClimaCore.Operators.SetValue(0.0))
gradc2f.(ψ) # fails

Please link to any related issues/PRs.

Describe the solution you'd like I think in the long term it would be good to allow fields as boundaries (assuming they are on the correct space) for all operators.

As noted, I think this is low priority.

kmdeck commented 2 years ago

@LenkaNovak @akshaysridhar may be interested in following; Lenka and I may take a stab at this after looking at Ben's PR where we added in spatially varying BC for divf2c also.