(Towards #868) Metadata support for fields with nlayers and ndata

[arporter]

No description provided.

[codecov]

Codecov Report

❌ Patch coverage is 65.21739% with 8 lines in your changes missing coverage. Please review.
✅ Project coverage is 99.94%. Comparing base (f0c5971) to head (892245f).

Files with missing lines	Patch %	Lines
src/psyclone/domain/lfric/lfric_arg_descriptor.py	61.11%	7 Missing ⚠️
src/psyclone/parse/kernel.py	80.00%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #3209      +/-   ##
==========================================
- Coverage   99.95%   99.94%   -0.02%     
==========================================
  Files         387      387              
  Lines       54320    54340      +20     
==========================================
+ Hits        54298    54310      +12     
- Misses         22       30       +8     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[arporter]

Please could @TeranIvy and/or @tommbendall (or someone else?) take a look at the updated documentation that this PR introduces. Once we're happy with that, I'll proceed to implementation.

[arporter]

Tagging @stevemullerworth so he can join the fun too :-)

[stevemullerworth]

One of my changes (for multi-data fields) suggests a change to the implementation. The others are document changes only, I think.

[stevemullerworth]

I would suggest avoiding the idea that there is a particular model mesh. For example, we will have apps that regrid data from one number of levels to another.

Suggested change

	Typically, a field will have the same number of vertical layers as the
	model mesh. However, this is not a requirement and the number of layers
	can be as few as one (a 2D field).
	While an application may model a set of fields that all have the same number of
	vertical levels, this is not a requirement. An application may have fields on a range
	of different levels including fields on a single level (a 2D field).

[stevemullerworth]

I may be misunderstanding this. I am assuming it relates to a kernel which, for example, takes an argument which may be, say, an 70-level field but could be a 50-level or 40-level field. But the kernel should only set 40 levels (or fewer) of the field. Most typically, perhaps, it would be useful when setting the bottom level of a field of any number of levels.

Suggested change

	that is not the same as the extruded mesh then this must be specified using
	that may be different from the mesh of the input field then this must be specified using

[arporter]

I don't know what the use case is I'm afraid. However, when we discussed with Chris and Iva we were definitely talking about fields that did not have e.g. 70 levels because that then changed the values in the dof map. On the other hand, if all fields actually have the same number of levels but only populate a subset of them then all of their dofmaps will be identical (for a given FS). I think this needs clarification from @TeranIvy or @tommbendall?

[arporter]

Coming back to this now, we recently agreed with @TeranIvy and @christophermaynard that the nlevels extension was the more straightforward and we would do this first. However, I've now re-discovered this thread and realised I still need an answer to the question: does a field with e.g. 40 levels actually have 70 but only makes use of the first 40? If so, it will have the same dofmap as the 70-level field.
Then there's the special case of 2D fields - perhaps this answers my first question as presumably you don't only make use of the 1st of 70 levels as that's a lot of wasted storage? Tagging @stevemullerworth and @tommbendall too.

[stevemullerworth]

Each field needs to have correct metadata in the kernel. The metadata can be used to identify the correct dofmap and the correct nlayers for the field. That's a strict rule.

An application can create fields with a 40-level extrusion and a 70-level extrusion (and a 2d field which in practice is a 1-level "extrusion"). Each extrusion requires a different 3D mesh and gives you a different nlayers, dofmap etc.

My suggested change refers to a kernel that, for example, applies fieldinput to fieldoutput. fieldoutput has 70 levels. fieldinput has 40 levels. So the kernel API could have two nlayers arguments and two sets of dofmap arguments. The two fields would be associated with different LFRic function spaces that are each attached to different 3D meshes and so different extrusions.

It's up to the kernel to know that it is looping over the fieldinput nlayers value and not the fieldoutput nlayers value.

Technically, one could write a kernel that loops over a subset of layers. But in that case, the number of layers would need to be a separate kernel input and the nlayers input would be ignored (except to check that we're not looping over more than nlayers levels).

[arporter]

Thanks Steve, I understand things a bit more now. Currently, we pass nlayers obtained from the first field/operator kernel argument:

I think what we're saying therefore is that any field/operator that has a number of levels that differs from the first field/operator argument must have nlayers=xx in its metadata. In turn, this will mean an additional dofmap is supplied to the kernel for each distinct combination of function space and nlayers. A kernel developer is of course free to specify nlayers for every field/operator kernel argument. That would just mean we would pass one of the dofmaps to the kernel twice but that's not a massive overhead (as opposed to requiring that all kernels be re-written).

[arporter]

Or, have I misunderstood and the proposal is to mandate that all kernel metadata specify nlayers for every field/operator argument?

[tommbendall]

Just catching up on this. A change to enforce nlayers metadata for every field/operator would be a very daunting task! So I would support the former proposal: only including an nlayers metadata argument for the fields/operators that differ from the first field/operator.

For existing kernels that have fields/operators on different meshes, we are currently managing this by declaring the fields/operators to be on different function spaces.

[arporter]

I now have:

If people are OK with that, the next question is whether to require nlayers to be an integer (in which case the 'labels' will have to be at least declared and initialised locally in the module containing the kernel) or whether we go for a string of some fixed length. That will avoid the need to declare anything but will require a character member of arg_type. As I've said elsewhere, I'm a bit uneasy about this because we've not done it before and there may be a good reason for that (and then again, it may just be because we've always wanted to restrict it to certain recognised values and integer parameters were an efficient way of enforcing that using Fortran alone).

[stevemullerworth]

For LFRic atmosphere, It would be useful to have the same flexibility for giving named values as is done for NLAYERS: in fact, it is more useful to have the flexibility here than for NLAYERS as there is less mixing and matching of layer numbers.

We have a lot of physics kernels with ~100 fields with a handful of multidata choices distinguished by ANY_DISCONTINUOUS_SPACE_1, ANY_DISCONTINUOUS_SPACE_2 etc. If they are converted to GH_RUNTIME we would need separate dofmaps etc. for each of the many fields rather than for each of the ANY_DISCONTINUOUS_SPACE choices. Example kernel here [MOSRS password protected]:

https://code.metoffice.gov.uk/trac/lfric_apps/browser/main/trunk/interfaces/physics_schemes_interface/source/kernel/bl_exp_kernel_mod.F90

[tommbendall]

Just commenting that I agree with this point! In fact I think the most common situation is that NDATA will be a runtime variable, but there might be many multidata fields being passed as arguments the same kernel, but all using the same NDATA value.

I'm not sure what the best solution to this is... because the NDATA values are science-related, I don't think we want any specific values hard-coded in arguments_mod.

Could we define a local variable for NDATA to take? e.g.

type(ndata_type) :: LAND_TILES
type(ndata_types) :: NUM_AEROSOLS

arg_type(GH_FIELD, GH_REAL, GH_READ, W2, NDATA=4)
type(arg_type) :: meta_args(4) = (/  &
     arg_type(GH_FIELD, GH_REAL, GH_READWRITE, W3),  &
     arg_type(GH_FIELD, GH_REAL, GH_READ,  W3, NDATA=LAND_TILES), &
     arg_type(GH_FIELD, GH_REAL, GH_READ,  WTHETA, NDATA=NUM_AEROSOLS), &
     arg_type(GH_FIELD, GH_REAL, GH_READ,  W3, NDATA=NUM_AEROSOLS)  &
/)

[stevemullerworth]

This is the sort of thing I had envisaged (though originally I was thinking of a direct but more meaningful replacement for ANY_*SPACE entries).

One of the challenges of moving away from ANY_DISCONTINUOUS_SPACE settings was deciding how and where to name scientifically-meaningful versions given that arguments_mod is in core. One could have a physics module to store them which would enable a common set to be used by several kernels.

I think they would need to be initialised so as not to cause any compiler warnings? If so, that's another reason for hiding them in a module.

[arporter]

Great, that's the same as I was imagining for NLAYERS so it's straightforward to do from a PSyclone point of view. I can see your issue with names in core and physics but I'll let you wrangle that :-)

[tommbendall]

Something else to consider that Tom M has just mentioned is that we have two types of multidata field. Those with ndata_first and those without. We may also need to capture that

[stevemullerworth]

Suggested change

	of layers which is not ``GH_RUNTIME`` then only one dofmap is passed to
	the kernel for those arguments.
	of layers which is not ``GH_RUNTIME`` then only one dofmap, the dofmap of the
	first field listed in the metadata, is passed to the kernel for those arguments.

[arporter]

I've had a play and the old metadata parser in PSyclone is quite happy to accept quoted strings, e.g.

arg_type(gh_field, gh_real, gh_read, w3, nlevels='double')

This means there's no need to have pre-declared and initialised integer values. However, argument_mod::arg_type in LFRic core would have to be updated to contain an nlevels:

 ...
     !> Optional metadata (fields only) for a type of stencil map to use.
     !! One of {XORY1D, X1D, Y1D, CROSS, REGION, CROSS2D}.
     integer :: stencil_map = -1
     !> Optional metadata (fields only) for inter-grid mapping kernels.
     !! One of {GH_FINE, GH_COARSE}.
     integer :: mesh_arg    = -1
     !> Optional metadata specifying the number of vertical levels for a field.
     character(len=SOME_LEN) :: nlevels=""
  end type arg_type

where SOME_LEN is some parameter with a suitable value. I'm not sure why this approach wasn't taken with other meta-data components so it could be there's some significant drawback to it?