Domain kernels and continuous fields

[iboutle]

I'm trying to pass the W2 wind into a DOMAIN kernel. PSyclone immediately tells me that this is not allowed, only discontinuous fields can be used in DOMAIN kernels.

I can see why this would cause issues if trying to write to a continuous field, however I'm passing it in as GH_READ, I only want to read from the field.

If I change the Psyclone metadata to define the function space as ANY_DISCONTINUOUS_SPACE instead of W2, it works and produces the results I would expect, i.e. the DOMAIN kernel can successfully read from the W2 function space.

Is PSyclone just being over-defensive in preventing a W2/GH_READ field being passed to a DOMAIN kernel, or is there good reason not to do this?

[Steve Mullerworth (stevemullerworth)]

From an LFRic developer point of view I can't think of reasons why this could not work as you have observed. However, I note some issues and questions below. The PSyclone developers would need to be consulted.

Can you summarise your requirement. Are you using both horizontal and vertical components of the W2 field?

Could you write this kernel as a non DOMAIN kernel - looping over cells? The DOMAIN kernels were originally provided with physics requirements in mind. A key aim was to allow the data order to be reordered from column-first to layer-first in line with the API of the physics schemes being transferred as is from the UM. If you need to reorder the genuine discontinuous fields in this kernel, can you split out the part of the operation that uses W2 into a separate kernel?

Does a DOMAIN kernel make sense with a W2 field - which traditionally holds an FEM field?

Note, pretending it is a discontinuous field may produce wrong answers if we ran without redundant computation (currently, we run with redundant computation as the reduction in halo swaps outweighs the cost of extra computation). In such a configuration, the W2 field on most ranks will contain annexed dofs that require a halo swap prior to the kernel call. Their marking as a discontinuous field may mean they don't get one. PSyclone would have to properly take this into account. Also, I note that a domain kernel declarations can reference the seg_len of input arrays as the PSy layer can be written to call the kernel in chunks with different segment sizes. The segment size of a W2 field would be different to the segment size of a W3 field and would need to be computed correctly.

[iboutle]

Thanks Steve,

I'm wanting to pass the horizontal wind components on their native grid into some of the existing physics code, hence the DOMAIN kernel.

A separate kernel would be possible but clunky - it would be doing very little actual work but require duplicating a fair amount of setup that is already happening in the physics routines I'm trying to pass it into.

Whether it makes sense possibly highlights an issue with our interchangeability at lowest-order between FE and FV fields - an FV field on the staggered c-grid certainly does make sense in a domain kernel, and W2 is the machinery by which we treat those fields throughout the model.

Your point about the annexed dofs is a good one. I seem to remember when trying to setup LFRic at BoM that there are plenty of other things that don't work if we turn off the redundant compute of the annexed dofs, but that's not a good reason to introduce more!

I suppose my other option is the broken W2 wind, which ought to work?!

[Mike Hobson (mike-hobson)]

I'm consious of not "teaching my grandmother to suck eggs", as I you've been around the LFRic code for a long time, but I'm intrigued about what you are doing here. You mention that you want the horizontal wind components on their native grid, but the wind field we store on its native grid doesn't contain "components" as such, but rather the fluxes across the faces that the dofs are sitting on. What we might call "x-" and "y-" components of our horizontal wind field are almost certainly not properly orthogonal or related to any geographical direction.

[iboutle]

Since you ask... we've recently discovered a nasty instability where we get quite strong zig-zagging (2 grid-point waves) in the horizontal wind components (in both horizontal and vertical directions). Unchecked these cause the model to fail.

The problem is that almost nothing actually sees the native grid winds - the transport and all of the physics operates on the W3 winds, where the averaging to obtain these removes the zig-zagging. Hence nothing can respond to control the noise.

We can exert some control over it by allowing the boundary layer scheme to see the vertical shear on the native grid, and respond to that. The orientation or non-orthogonality doesn't really matter here, all I'm interested in is getting some estimate of the vertical shear on the native grid to enable the vertical mixing to act upon and limit/remove it.

[Mike Hobson (mike-hobson)]

Thanks for that. Sounds like fun!

[Steve Mullerworth (stevemullerworth)]

I feel it would be preferable to move the components of the wind that you need onto a W2h broken space to preserve the behaviour when running without redundant computation and to avoid the need to stretch the requirement for domain kernels too far.

For interest, are you transposing both W2 and cell-centre fields to layer-first ordering before passing to the physics?

Your comment does highlight a need to support build-time configuration of PSyclone so that we can test with and without redundant computation. It sounds like we are building technical debt into the model that would hurt us if for any reason we did not want to do redundant computation (though maybe that is unlikely...?).

[iboutle]

I've not got as far as reordering it yet - I don't think it's problematic to do, but am still pondering the best approach given the comments above.

I can't imagine it will ever be computationally preferable to halo-swap the annexed dofs rather than compute them, but as you say, that doesn't mean it shouldn't work!