Climatemachine.jl: Question: Mapping equations to functions

[charleskawczynski]

This is more of a question / clarification than an issue

I made a gist of a compressed/pseudo code of DGModel_kernels, and method calls within these kernels to functions defined in the AtmosModel. I would like to review/clarify how these methods are expressing the equations since they're used in multiple places:

Let G = ∇Y

∂Y
-- = - ∇ • F_{nondiff}(Y) - ∇ • F_{diff}(Y,G) + S(Y)
∂t

| Equation term / purpose | Method |
|---------------------------------------:|:----------------------------------------|
| Pre-compute auxstate | update_aux! |
| Pre-compute gradients | gradvariables! |
| ? | volumeviscterms! -> diffusive! |
| ? | faceviscterms! -> diffusive! |
| F_{nondiff}(Y) | volumerhs! -> flux_nondiffusive! |
| F_{diff}(Y,G) | volumerhs! -> flux_diffusive! |
| S(Y) | volumerhs! -> source! |

Concretely, I'd like to fix #701 so that subsidence is a source, not a flux. However, the expression for subsidence involves a gradient, so it seems that source! args need to include state::Grad. I'd like to understand how the other methods called enter this equation, or what their purpose is.

[charleskawczynski]

cc: @lcw @jkozdon should I add a loop to compute gradients prior to the source! call to fix this? Then I can update the table to say S(Y) -> S(Y,G), correct?

[jkozdon]

You should need to add another loop to compute the gradient, it's already computed and stored in what you call G. You just need to pass this to the source.

That said, if you shove too much stuff into source via this mechanism you could end up with some numerical stability problems. For instance, you wouldn't want to put terms that should really be in F_{nondiff}(Y) into S.

[jkozdon]

| Equation term / purpose | Method |
|---------------------------------------:|:----------------------------------------|
| Pre-compute auxstate | update_aux! |
| transform from state to variables to take gradient of (typically primitive) | gradvariables! |
| linear transform of gradient variables (for instance to store 6 stresses instead of 9 gradients variables). These would be the DG auxiliary variables | faceviscterms! / volumeviscterms! -> diffusive! |
| F_{nondiff}(Y) | volumerhs! -> flux_nondiffusive! |
| F_{diff}(Y,G) | volumerhs! -> flux_diffusive! |
| S(Y) | volumerhs! -> source! |

[jkozdon]

That said, if you shove too much stuff into source via this mechanism you could end up with some numerical stability problems

If you do have problems, if you can give some minimal PDE (something akin to convection-diffusion) that has similar behavior we can see how to implement things in a (more) stable manner.

[charleskawczynski]

Thanks! This is super helpful.

You should need to add another loop to compute the gradient

Do you mean _shouldn't_?

That said, if you shove too much stuff into source via this mechanism you could end up with some numerical stability problems. For instance, you wouldn't want to put terms that should really be in F_{nondiff}(Y) into S.

Good to note. Is there a general rule on the type of sources that would result in issues? (i.e., does this mostly pertain to stiff sources?). Regardless, it sounds like we can alleviate thinking about what might be problematic by simply re-formulating sources with gradients to flux terms.

The particular source term is of the form

∂Y
--- = - ∇ • F_{nondiff}(Y) - ∇ • F_{diff}(Y,G) - W ẑ ⋅ ∇Y
∂t

Where W = -Dz

[jkozdon]

Do you mean _shouldn't_?

Yes

Good to note. Is there a general rule on the type of sources that would result in issues? (i.e., does this mostly pertain to stiff sources?)

I don't know because I have never tried to put a term involving the derivative in the source.

If it should be just non-conservative form that's different than moving it to the source in general. @fxgiraldo might have more experience with non-conservative forms, but there is some DG work related to this that could be looked at.

[simonbyrne]

To address Charlie's question at the start, I haven't quite thought this through, but one thing I have been kicking around in my head as part of #690 is if we could define "DG gradient" and "DG divergence" operators so that we could write the main DG loop as something like:

@fuse begin
  transform = gradvariables.(state, aux, t)
  ∇transform = dg_gradient(transform)
  diffstate = diffusive.(∇transform, state, aux, t)
end

@fuse begin
  F_nondiff = flux_nondiff.(state, aux, t)
  F_diff = flux_diff.(state, diffstate, aux, t)
  S = source.(state, aux, t)
  rhs = dg_divergence(F_nondiff) .+ dg_divergence(F_diff) .+ S
end

where state, aux etc are "Field" objects i mentioned in #690 (obviously still need a lot more thought on how boundary conditions and numerical fluxes are handled, and how @fuse would work).

Then users could chop and change what they need. e.g. to deal with the needs of the ocean model, they could have an extra 2D field, which would be handled automatically via broadcasting. update_aux! would be replaced by additional code between @fuse blocks.

[blallen]

@simonbyrne, @sandreza has wanted to do this forever, e.g. we think this is a great idea and makes it easier to dispatch on numerical flux and stuff

[tapios]

Just to clarify a few issues raised above (@jkozdon, @simonbyrne):

• The term to be added is a nonconservative source; it cannot be written as a flux divergence.
• The source does not pose stiffness problems (e.g., W is fixed and smooth; the gradients involved are slowly varying); in fact, we have computed it in the way we want it here in LES for decades.
• The only alternative to the straightforward imposition of a source -W ẑ . gradient(q) is to add -divergence(W ẑ q) + q. divergence(W ẑ), with the second term becoming a source. This would work in DYCOMS because divergence(W ẑ) is constant and can easily be precomputed. But it would not work in more general circumstances where we need this functionality.

[jkozdon]

The source does not pose stiffness problems (e.g., w_sub is fixed and smooth; the gradients involved are slowly varying); in fact, we have computed it in the way we want it here in LES for decades.

When I said numerical instabilities before I wasn't thinking stiffness but that a positive real eigenvalue could be introduced in the ODE system being solved.

To clarify, I don't know that this will happen but just that if weird instabilities start popping up when subsidence is on that this could be the culprit. (It might indeed have no problems, I don't know.)

[tapios]

The source does not pose stiffness problems (e.g., w_sub is fixed and smooth; the gradients involved are slowly varying); in fact, we have computed it in the way we want it here in LES for decades.

When I said numerical instabilities before I wasn't thinking stiffness but that a positive real eigenvalue could be introduced in the ODE system being solved.

To clarify, I don't know that this will happen but just that if weird instabilities start popping up when subsidence is on that this could be the culprit. (It might indeed have no problems, I don't know.)

The physics here is that subsidence is stabilizing. (Subsidence is what ultimately balances surface fluxes and turbulent entrainment in the boundary layer mass and energy budgets.)