CliMA / CliMA/ClimaCore.jl

Cannot use DivergenceC2F with SetValue boundary conditions in an operator matrix

Open
#2,640 0 comments 0 reactions 0 assignees View on GitHub
bug
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

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.