Cannot use DivergenceC2F with SetValue boundary conditions in an operator matrix
- Dominant language
- Julia
- Stars
- 117
- Forks
- 19
- Avg merge
- 3d 4h
- Merged PRs (30d)
- 41
Description
## Describe the bug
Using the `DirichletC2F` operator with a `SetValue` boundary condition inside implicit timesteppers seems to be no longer possible, since the call to `operator_matrix()` fails.
I believe this bug is a result of the FD operator revamp, #2544, and thus present since v0.16.0.
The [docs section on boundary conditions](https://clima.github.io/ClimaCore.jl/stable/howto/boundary_conditions/#Vertical-boundaries-(finite-difference-operators)) explicitly mentions that `SetValue` boundary conditions are supported for `DirichletC2F`, so this is unexpected.
_Side note:_ I have tried using an equivalent `SetGradient` condition on the `GradientF2C` operator instead.
This silently leads to `NaN` values in the operator matrix and the solution becomes `NaN` when using an implicit timestepping method.
I believe this change of behavior can also be inferred from the [release notes for 0.16.0](https://github.com/CliMA/ClimaCore.jl/blob/main/NEWS.md#v0160).
## Reproducer
I have adapted this from [the cubed sphere example from the docs](https://clima.github.io/ClimaCore.jl/stable/tutorials/extruded_sphere/#3.-Explicit-horizontal-and-implicit-vertical-tendencies) (by the way, very nice with the new structure!):
```julia
using ClimaComms
ClimaComms.@import_required_backends
import ClimaCore: Domains, Meshes, Spaces, Fields, Geometry, Operators
import ClimaCore.MatrixFields
import ClimaCore.MatrixFields: @name, ⋅
vert_domain = Domains.IntervalDomain(
Geometry.ZPoint(0.0),
Geometry.ZPoint(1.0);
boundary_names = (:bottom, :top),
)
vert_mesh = Meshes.IntervalMesh(vert_domain; nelems = 10)
device = ClimaComms.device()
vert_space = Spaces.FaceFiniteDifferenceSpace(device, vert_mesh)
C3 = Geometry.WVector
divᵥ = Operators.DivergenceC2F(;
bottom = Operators.SetValue(C3(0.0)),
top = Operators.SetValue(C3(0.0)),
)
divᵥ_matrix = MatrixFields.operator_matrix(divᵥ)
```
Error message (latest version of the main branch, corresponding to v0.16.2):
```julia
ERROR: MethodError: no method matching operator_matrix(::ClimaCore.Operators.DirichletOperator{ClimaCore.Operators.DivergenceC2F, @NamedTuple{…}})
The function `operator_matrix` exists, but no method is defined for this combination of argument types.
```
Contributor guide
No contributing guide indexed for this repository
Research direction
Start by running the reproducer through MatrixFields.operator_matrix, then inspect the DirichletC2F and SetValue path for operator matrices. Compare the behavior with the boundary-conditions documentation and the v0.16.0 release notes. Done means the documented SetValue case no longer raises MethodError and the related operator-matrix path does not produce NaN values.
Written by the indexing model from the issue text.
Assessment
- Tech stack
- julia
- Domain
- backend
- Issue type
- Bug
- Difficulty
- 3/5
- Estimated time
- 1-2 days
- Activity status
- Active
- Clarity
- Mostly clear
- Newbie friendliness
- 68/100