Spatially varying BC dont work for grads?

[kmdeck]

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

Then, let's try and apply value as a boundary flux to divergence:

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

Another scenario is a spatially varying state BC - i.e. apply it as a boundary to the gradient:

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

Any idea why spatially varying BC doent seem to work for gradients?

[kmdeck]

@charleskawczynski @juliasloan25

Charlie, more generally, it might be good to chat about how best to apply dirichlet conditions. We need them for tests. Right now, we want to convert a (field of) state BC values to a flux, using finite difference: flux_bc \approx -(state_bc - state_center)/dz

Then the RHS code always handles boundary fluxes, and applies boundary conditions in the divergence operator.

This is creating lots of kind of hacky code, because the state_bc and state_center may be on different spaces, the flux_bc needs to be on the face space at the surface, etc.

Im wondering if we should be making use of applying spatially varying BC to the gradient directly instead (dispatching off of the type of BC to create the operators, i.e. with BC on gradc2f or on divf2c depending on if it is a state or flux bc).

[charleskawczynski]

Sorry for the late response, it does seem that doing this bycolumn works, though!

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

function get_center_space(::Type{FT}) where {FT}
    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 = Domains.IntervalDomain(
        Geometry.XPoint(xlim[1]),
        Geometry.XPoint(xlim[2]);
        periodic = true,
    )
    domain_y = Domains.IntervalDomain(
        Geometry.YPoint(ylim[1]),
        Geometry.YPoint(ylim[2]);
        periodic = true,
    )
    plane = Domains.RectangleDomain(domain_x, domain_y)

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

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

    return Spaces.ExtrudedFiniteDifferenceSpace(horzspace, vert_center_space)
end

FT = Float64
center_space = get_center_space(FT)

ψ = Fields.ones(center_space)

surface_field = Fields.zeros(axes(Spaces.level(ψ, 1)))
value = surface_field # same type/instance of underlying space as horizontal space of \psi

gradc2f_no_bc = Operators.GradientC2F()
divf2c = Operators.DivergenceF2C(
    top = Operators.SetValue(Geometry.WVector.(value)),
    bottom = Operators.SetValue(Geometry.WVector(FT(0))),
)
@. divf2c(gradc2f_no_bc(ψ)) # runs

# gradc2f = Operators.GradientC2F(
#     top = Operators.SetValue(value),
#     bottom = Operators.SetValue(FT(0))
# )
# gradc2f.(ψ) # fails

Fields.bycolumn(axes(ψ)) do colidx
    value_col = value[colidx]
    gradc2f_col = Operators.GradientC2F(
        top = Operators.SetValue(value_col),
        bottom = Operators.SetValue(FT(0))
    )
    gradc2f_col.(ψ) # works
end

[kmdeck]

thank you, Charlie! Im hoping switching to column wise operations makes the code a little cleaner.