Update vector hyperdiffusion to act on full vectors

[bischtob]

Currently the vector hyperdiffusion operation in ClimaAtmos does not operate like the standard vector hyperdiffusion. This is because the current vector Laplacian used in the implementation does not act on three-dimensional vector fields.

Loosely speaking, we are currently doing this:

scurl = Operators.Curl()
sdiv = Operators.Divergence()
wcurl = Operators.WeakCurl()
wgrad = Operators.WeakGradient()

# typeof(u_h) -> Covariant12Vector
laplacian(u_h) = @. wgrad(sdiv(u_h)) - Covariant12Vector(wcurl(Covariant3Vector(scurl(u_h))))
# typeof(laplacian(u_h)) -> Covariant12Vector

We should use the correct vector Laplacian by chaining together differential operators and coordinate transforms that act on three-dimensional vector fields. For example like this:

scurl = Operators.Curl()
sdiv = Operators.Divergence()
wcurl = Operators.WeakCurl()
wgrad = Operators.WeakGradient()

# typeof(u_bar) -> Covariant123Vector
laplacian(u_bar) = @. wgrad(sdiv(u_bar)) - Covariant123Vector(wcurl(Covariant123Vector(scurl(u_bar))))
# typeof(laplacian(u_bar)) -> Covariant123Vector

# TODO: 
Interpolate third component to faces.

In this pseudocode, the main software issue lies with scurl(u_bar) because our spectral element curl operator currently does not operate on Covariant123Vector types out of the box. The fix seems reasonably simple in ClimaCore.jl, but a hacky solution also exists.

[bischtob]

@tapios, let me know if this is how you think we should address the vector hyperdiffusion issue or if I am misunderstanding the problem.

It seems that after the operation is done, we just need to extract the correct components for the tendencies.

[bischtob]

cc. @akshaysridhar

[tapios]

@tapios, let me know if this is how you think we should address the vector hyperdiffusion issue or if I am misunderstanding the problem.

It seems that after the operation is done, we just need to extract the correct components for the tendencies.

We need some more detail here:

1. The derivatives should be constrained to be horizontal (we do not want to take derivatives across horizontal coordinate surfaces--think of a highly anisotropic diffusion tensor sandwiched in here that has no component across horizontal coordinate surfaces).

2. We need the details on how the vertical interpolations are done, using u_bar following https://clima.github.io/ClimaAtmos.jl/dev/equations/#Auxiliary-and-derived-quantities

[bischtob]

Thank you for the clarifications.

We need some more detail here:

1. The derivatives should be constrained to be horizontal (we do not want to take derivatives across horizontal coordinate surfaces--think of a highly anisotropic diffusion tensor sandwiched in here that has no component across horizontal coordinate surfaces).

In the code, we use the spectral element differential operators only for hyperdiffusion (e.g., there should only be derivatives wrt the first two generalized coordinates here).

2. We need the details on how the vertical interpolations are done, using u_bar following https://clima.github.io/ClimaAtmos.jl/dev/equations/#Auxiliary-and-derived-quantities

100%. We will diffuse u_bar (e.g., a Covariant123Vector if I am not mistaken).

The output from this operation would be a Covariant12Vector for the divergence term and a Covariant123Vector for the curl term, aka the second term would have a component in the direction of e^3, e.g. it is the second term that is producing a tendency for u_3.

Let me know if this matches your idea.

[OsKnoth]

As I understand correctly we want to add hyperdiffusion for the w (u_3) component. We can do this by seeing w as a scalar and have to apply the machinery for scalar hyperdiffusion to w with the right metric terms (have to be interpolated from cell to faces)

[bischtob]

@OsKnoth, I guess that is possible in theory, but I think @tapios suggested to actually do a vector Laplacian version of hyperdiffusion. Seems more elegant that way to me.

[bischtob]

Technical comment:

ClimaCore.jl does not provide a horizontal curl (spectral element operator) that acts on Covariant123Vectors. A simple, workaround in pseudocode may look like this:

# typeof(u_bar) -> Covariant123Vector
omega = scurl(Covariant12Vector(u_bar)) + scurl(Covariant3Vector(u_bar))
curlcurl(u_bar) = wcurl(Covariant12Vector(omega)) + wcurl(Covariant3Vector(omega)) # -> Contravariant123Vector

cc. @akshaysridhar

[akshaysridhar]

Technical comment:

ClimaCore.jl does not provide a horizontal curl (spectral element operator) that acts on Covariant123Vectors. A simple, workaround in pseudocode may look like this:

# typeof(u_bar) -> Covariant123Vector
omega = scurl(Covariant12Vector(u_bar)) + scurl(Covariant3Vector(u_bar))
curlcurl(u_bar) = wcurl(Covariant12Vector(omega)) + wcurl(Covariant3Vector(omega)) # -> Contravariant123Vector

cc. @akshaysridhar

Expanding this following REPL tests with existing ClimaCore Operators - Field values are truncated in the MWE below (given u_bar is accessible from cache variables): χ denote intermediate variables in this demonstration.

julia> u_bar = p.ᶜu
ClimaCore.Geometry.Covariant123Vector{Float32}-valued Field
julia> u_bar_12 = @. Geometry.Covariant12Vector(u_bar)
ClimaCore.Geometry.Covariant12Vector{Float32}-valued Field:
julia> u_bar_3 = @. Geometry.Covariant3Vector(u_bar)
ClimaCore.Geometry.Covariant3Vector{Float32}-valued Field:
julia> χ₁ = @. Geometry.Contravariant123Vector(scurl(u_bar_12)) + Geometry.Contravariant123Vector(scurl(u_bar_3))
julia> χc₁₂ = @. Geometry.Covariant12Vector(χ₁)
julia> χc₃ = @. Geometry.Covariant3Vector(χ₁)
julia> result_contra123 = @. Geometry.Contravariant123Vector(wcurl(χc₁₂)) + Geometry.Contravariant123Vector(wcurl(χc₃))
julia> result_cov123  = @. Geometry.Covariant123Vector(result_contra123)

i.e. the ClimaCore block contains the 2 operations we need separately (as in the "hack" above), an additional elseif RT <: Contravariant123Vector block will complete the curl operator specification for the Covariant123Vector type.

[bischtob]

Ok. I will open a draft PR that builds out the "hack" first. We can spec out how we would replace the hack with a new ClimaCore operations and open a separate PR there. This will allow us to run out test cases and spec out new ClimaCore functionality at the same time.

[akshaysridhar]

Ok. I will open a draft PR that builds out the "hack" first. We can spec out how we would replace the hack with a new ClimaCore operation and open a separate PR there. This will allow us to run out test cases and spec out new ClimaCore functionality at the same time.

Noted, thanks!

[bischtob]

We will also compare the math to the case of treating u_3 as a scalar wrt hyperdiffusion per @OsKnoth's point.

[tapios]

Technical comment:
ClimaCore.jl does not provide a horizontal curl (spectral element operator) that acts on Covariant123Vectors. A simple, workaround in pseudocode may look like this:

# typeof(u_bar) -> Covariant123Vector
omega = scurl(Covariant12Vector(u_bar)) + scurl(Covariant3Vector(u_bar))
curlcurl(u_bar) = wcurl(Covariant12Vector(omega)) + wcurl(Covariant3Vector(omega)) # -> Contravariant123Vector

cc. @akshaysridhar

Expanding this following REPL tests with existing ClimaCore Operators - Field values are truncated in the MWE below (given u_bar is accessible from cache variables): χ denote intermediate variables in this demonstration.

julia> u_bar = p.ᶜu
ClimaCore.Geometry.Covariant123Vector{Float32}-valued Field
julia> u_bar_12 = @. Geometry.Covariant12Vector(u_bar)
ClimaCore.Geometry.Covariant12Vector{Float32}-valued Field:
julia> u_bar_3 = @. Geometry.Covariant3Vector(u_bar)
ClimaCore.Geometry.Covariant3Vector{Float32}-valued Field:
julia> χ₁ = @. Geometry.Contravariant123Vector(scurl(u_bar_12)) + Geometry.Contravariant123Vector(scurl(u_bar_3))
julia> χc₁₂ = @. Geometry.Covariant12Vector(χ₁)
julia> χc₃ = @. Geometry.Covariant3Vector(χ₁)
julia> result_contra123 = @. Geometry.Contravariant123Vector(wcurl(χc₁₂)) + Geometry.Contravariant123Vector(wcurl(χc₃))
julia> result_cov123  = @. Geometry.Covariant123Vector(result_contra123)

i.e. the ClimaCore block contains the 2 operations we need separately (as in the "hack" above), an additional elseif RT <: Contravariant123Vector block will complete the curl operator specification for the Covariant123Vector type.

Is this only involving horizontal derivatives? It looks to me from https://clima.github.io/ClimaCore.jl/stable/operators/#ClimaCore.Operators.Curl that scurl also has the vertical derivatives, which we do not want here.

[bischtob]

Is this only involving horizontal derivatives? It looks to me from https://clima.github.io/ClimaCore.jl/stable/operators/#ClimaCore.Operators.Curl that scurl also has the vertical derivatives, which we do not want here.

I believe the docs are not correct, because in the code here https://github.com/CliMA/ClimaCore.jl/blob/77b6294fd9f55dba925b8b28bb8c5ca49de6dd68/src/Operators/spectralelement.jl#L1108 there seem to only be spectral derivatives in the 1 & 2 directions.

[tapios]

Is this only involving horizontal derivatives? It looks to me from https://clima.github.io/ClimaCore.jl/stable/operators/#ClimaCore.Operators.Curl that scurl also has the vertical derivatives, which we do not want here.

I believe the docs are not correct, because in the code here https://github.com/CliMA/ClimaCore.jl/blob/77b6294fd9f55dba925b8b28bb8c5ca49de6dd68/src/Operators/spectralelement.jl#L1108 there seem to only be spectral derivative in the 1 & 2 directions.

Ok. It would be good opportunity then to make sure that operator names more clearly reflect what they do (e.g., we have some with h subscript for horizontal).

@simonbyrne Could you please review plans here and amend where needed, so we are using the correct operators (with horizontal derivatives)?

[bischtob]

Technical comment:
There seems to also be an issue with DSS not accepting Covariant123Vector fields as input.

ERROR: LoadError: AssertionError: S <: supportedvectortypes