Clarify default boundary conditions for `GradientC2F`

[valentinaschueller]

In §2.2 of the tutorial in the docs, the following is stated:

Center to face gradient

Uses the same stencil, but doesn't work directly:

sinz = sin.(column_center_coords.z)
gradc2f = ClimaCore.Operators.GradientC2F()
# ∇sinz = gradc2f.(sinz) ## this would throw an error

This throws an error because face values at the boundary are not well-defined:

I am on v0.14.50, and this does not seem to be true, ∇sinz = gradc2f.(sinz) will not throw an error.
I have tried this with some experiments and it seems that the following boundary condition is applied:
At indices 0, n+1, homogeneous Dirichlet data is assumed, that is, at (nonexisting) cell centers outside the domain.¹
This leads to a boundary gradient value that is half the size of what you would get by applying homogeneous Dirichlet boundary conditions the standard way (i.e., SetValue(0.0)).
If I understand the code behavior correctly, it seems that with SetValue(), the Dirichlet data is placed at the faces (indices 1/2, n+1/2).
Below I have included a script that I wrote to understand what's going on, maybe it is helpful.

I am not so familiar with the code, so I am not sure that this is actually what's going on.
But perhaps you can clarify this in the documentation/the tutorial.

---

Example script:

import ClimaCore as CC
using ClimaComms
using ClimaCorePlots
using Plots

function get_coord(lower_boundary::Float64, upper_boundary::Float64, nelems::Int)
    device = ClimaComms.device()
    domain = CC.Domains.IntervalDomain(
        CC.Geometry.ZPoint(lower_boundary),
        CC.Geometry.ZPoint(upper_boundary);
        boundary_names=(:bottom, :top),
    )
    mesh = CC.Meshes.IntervalMesh(domain, nelems=nelems)
    cspace = CC.Spaces.CenterFiniteDifferenceSpace(device, mesh)
    coord = CC.Fields.coordinate_field(cspace)
    return coord
end

function plot_gradient(lower_boundary::Float64, upper_boundary::Float64, nelems::Int, color)
    coord = get_coord(lower_boundary, upper_boundary, nelems)
    sinz = sin.(coord.z)

    bottom_bc = CC.Operators.SetValue(0.0)
    top_bc = CC.Operators.SetValue(0.0)
    gradc2f = CC.Operators.GradientC2F(bottom=bottom_bc, top=top_bc)
    # uncomment the line below to switch to GradientC2F without boundary conditions
    # gradc2f = CC.Operators.GradientC2F()
    ∇sinz = gradc2f.(sinz)
    @info parent(CC.Geometry.WVector.(∇sinz))[end]
    plot!(
        map(x -> x.w, CC.Geometry.WVector.(∇sinz)), 
        color=color, 
        label="nelems = $nelems",
        xlabel="z",
        ylabel="∇sin(z) approximation"
    )
end

lower_boundary = 0.0
upper_boundary = 10.0
nelems = 10

plot()
plot_gradient(0.0, 10.0, 10, :black)
plot_gradient(0.0, 10.0, 5, :blue)
plot_gradient(0.0, 10.0, 20, :green)
plot_gradient(0.0, 10.0, 100, :red)
savefig("grad_comparison.png")

Footnotes

1. I am following the numbering convention from here ↩

[imreddyTeja]

@valentinaschueller I agree that this is confusing and needs to be better documented. When I reproduced this with your provided example script, I also did not get an error. However, when I reran the script with bounds checking enabled, I did get an error. To enable bounds checks, start Julia with --check-bounds=yes.

When SetValue() is provided for both bottom and top, the bottom and top faces compute the gradient between themselves and the nearest interior center. The bottom and top faces take their values from the SetValue.

When no boundary conditions are applied, the top and bottom faces are treated the same as the interior faces. This means the bottom face tries to access the 0th index, and the top tries to access the n+1 index. In the script you sent, it is likely that there is nothing immediately adjacent to sinz in memory, so when the top and bottom face access invalid memory, the values are likely close to zero. Here is a slight modification to your script that shows putting something adjacent to sinz in memory, which gives different results

function plot_gradient(lower_boundary::Float64, upper_boundary::Float64, nelems::Int, color)
    coord = get_coord(lower_boundary, upper_boundary, nelems)
    # change `a` and `c`, and the result will change if no bcs are provided
    bigger_field = fill((a=-2.0, b=0.0, c=2.0), axes(coord.z))
    sinz = bigger_field.b
    @. sinz = sin(coord.z)
    bottom_bc = CC.Operators.SetValue(0.0)
    top_bc = CC.Operators.SetValue(0.0)
    gradc2f = CC.Operators.GradientC2F(bottom=bottom_bc, top=top_bc)
    # uncomment the line below to switch to GradientC2F without boundary conditions
    gradc2f = CC.Operators.GradientC2F()
    ∇sinz = gradc2f.(sinz)
    plot!(
        map(x -> x.w, CC.Geometry.WVector.(∇sinz)),
        color=color,
        label="nelems = $nelems",
        xlabel="∇sin(z) approximation",
        ylabel="z"
    )
end

[valentinaschueller]

Thanks @imreddyTeja for the explanation, very interesting!

From a user standpoint, I think I would appreciate a meaningful error if I missed to provide a BC.

But in case that is nontrivial to implement or affects performance, I think it makes sense to mention this in the docs.
Maybe something like:

If no boundary condition is provided, boundary faces will be treated the same way as interior faces.
This can lead to undefined behavior (if bounds checking is enabled, it will throw a BoundsError).

in the dev docs.
In the tutorial, I would suggest removing the note that it will throw an error.