Psyclone creating incorrect ordering of dofmaps in kernel argument list?

[MichaelWhitall]

Hi all,

I have a branch which (among other things) is adding a new kernel to perform sanity-checks on the prognostic precip fraction used with CoMorph:
vn3.1_max_in_cloud_checks
(see this issue page for more details: MetOffice/lfric_apps#243).

Currently if I compile and run with this branch I get failures due to referencing exner out-of-bounds in my new kernel lsp_precfrac_checks_kernel_mod.F90, which is invoked from cld_alg_mod.x90.

After a load of print statements, and hacking rose-stem so that it keeps the pre-processed source output by psyclone, I've finally found that the bug resides in some of the autogenerated code :( In my kernel's argument list I have:

subroutine lsp_precfrac_checks_code( nlayers,                                  &
                                     exner_w3,                                 &
                                     mv_wth, ml_wth, mi_wth,                   &
                                     ms_wth, mr_wth, mg_wth,                   &
                                     precfrac_wth,                             &
                                     ndf_wth, undf_wth, map_wth,               &
                                     ndf_w3,  undf_w3,  map_w3 )

Where this is to be called from cld_alg, I have:

      ! Sanity-checks on the prognostic precip fraction
      call invoke( lsp_precfrac_checks_kernel_type( exner,                     &
                                                    mr(imr_v),                 &
                                                    mr(imr_cl),                &
                                                    mr(imr_ci),                &
                                                    mr(imr_s),                 &
                                                    mr(imr_r),                 &
                                                    mr(imr_g),                 &
                                                    precfrac ) )

(note the arguments ndf_wth, undf_wth, map_wth, ndf_w3, undf_w3, map_w3 are intentionally missing as psyclone is meant to put them in automatically when it converts the above into a call to the subroutine lsp_precfrac_checks_code).

But if I look at the autogenerated code made by psyclone which actually calls this subroutine (the file called algorithm/cld_alg_mod_psy.f90), what it has autogenerated is:

    ! Call kernels and communication routines
    !$omp parallel default(shared) private(cell)
    !$omp do schedule(static)
    do cell = loop0_start, loop0_stop, 1
      call lsp_precfrac_checks_code(nlayers_exner, exner_data, mr_data, mr_1_data, mr_2_data, mr_3_data, mr_4_data, mr_5_data, &
&precfrac_data, ndf_w3, undf_w3, map_w3(:,cell), ndf_wtheta, undf_wtheta, map_wtheta(:,cell))
    enddo
    !$omp end do
    !$omp end parallel

Comparing this with my kernel argument list above, psyclone has wrongly changed the ordering of the arguments. I have ndf_w3, undf_w3, map_w3 at the end of the argument list in subroutine lsp_precfrac_checks_code, but in the call from the autogenerated code these arguments have been moved to before ndf_wtheta. Hence inside my subroutine, the dofmap which I've labelled as map_w3 is actually pointing at map_wtheta, which has one extra level in each column, causing my exner array to be vertically misaligned and be referenced out-of-bounds.

Anyone have any idea what's going on / how I can fix this? Is it a bug in psyclone? Or does psyclone just have a pre-defined order in-which the dofmaps must appear in the argument list, and I need to change my kernel interface to have the _w3 dofmap first, to match what it expects?? (some kernels use quite a few different dofmaps; if they must appear in a specific order in the argument list, what is that order / where is this documented??)

Cheers!
Mike

[Andrew Coughtrie (andrewcoughtrie)]

The order of the ndf, undf arguments is defined by metadata at the top of the kernel. So in your case the W3 field comes and first and then Wtheta so the ndf_w3 etc. will come before ndf_wth. The kernel section of the documentation describes this. There is also the kernel stub generator (the doc linked calls it genkernelstub, this is incorrect and it is now called psyclone-kern) which can be used to create the subroutine argument list for a kernel from the metadata.

[MichaelWhitall]

Hi thanks Andrew;

Re the metadata at the top of the kernel, do you mean this thing?:

type, public, extends(kernel_type) :: lsp_precfrac_checks_kernel_type
  private
  type(arg_type) :: meta_args(8) = (/                                          &
       arg_type(GH_FIELD, GH_REAL, GH_READ,  W3),     & ! exner_w3
       arg_type(GH_FIELD, GH_REAL, GH_READ,  WTHETA), & ! mv_wth
       arg_type(GH_FIELD, GH_REAL, GH_READ,  WTHETA), & ! ml_wth
       arg_type(GH_FIELD, GH_REAL, GH_READ,  WTHETA), & ! mi_wth
       arg_type(GH_FIELD, GH_REAL, GH_READ,  WTHETA), & ! ms_wth
       arg_type(GH_FIELD, GH_REAL, GH_READ,  WTHETA), & ! mr_wth
       arg_type(GH_FIELD, GH_REAL, GH_READ,  WTHETA), & ! mg_wth
       arg_type(GH_FIELD, GH_REAL, GH_WRITE, WTHETA)  & ! precfrac_wth
       /)
   integer :: operates_on = CELL_COLUMN
contains
  procedure, nopass :: lsp_precfrac_checks_code
end type

But the dofmap stuff doesn't appear here; only the array arguments. Do you mean the order in-which the dofmaps must appear is set by the ordering of the fields that are defined on those dofmaps? So in my case, the _w3 dofmap must appear first because the first array argument is defined in w3 space, whereas the _wtheta dofmap must come next since the next array is defined in wtheta space?

It strikes me that it'd be really easy to get this wrong without noticing that its wrong. e.g. if I added a new argument to a kernel with an already-long argument list using many different dofmaps, I might inadvertently change the order in-which each dofmap's corresponding space first appears in the list of input arrays. This will sometimes cause the model to fail, but it could potentially just scramble some of the data going in / out in a non-fatal way, so the model would perform worse with no way of telling why. Two things jump to mind:
a) This is exactly the sort of bug that you can catch by testing bit-comparison on different processor decompositions. Why was this key safety-net quietly dropped in the move from um to lfric? We urgently need to implement a bit-comp test for the physics in lfric_apps rose-stem!
b) Wouldn't it be easier / safer if the appearance of the dofmaps in a kernel's _code subroutine argument list was autogenerated by psyclone, consistent with the autogenerated call to that subroutine from the algorithm? Then the possibility of mere mortal human code-developers like me getting it wrong would be removed?

Cheers!
Mike

[Andrew Coughtrie (andrewcoughtrie)]

Yes that's the metadata (sorry if that wasn't as clear as it could be).

It can be tricky to get things right yes and we are always open to suggestions on improving the checks and tests that ensure issues don't creep in.

In answer to b), if I understand you correctly this is what the stub generator does, it takes the metadata and creates the kernels _code subroutine argument list and declarations for you. The metadata is what is used to tell PSyclone what it is going to be given through the invoke and what it should then extract from various places and send through to the kernel (e.g. stencils). For more information on how the LFRic DSL works see here: https://psyclone.readthedocs.io/en/stable/user_guide/lfric.html there is quite a lot of complexity in what can be done with PSyclone and Kernels.

On an additional note: it is standard practice to put the out/inout (GH_WRITE/GH_READWRITE/GH_INC etc.) fields first in the kernel metadata, this would result in the Wtheta arguments coming before the W3 ones in your case as your only WRITE field is on Wtheta. Info on the codeing standards for writing kernels is here.

[MichaelWhitall]

PS thanks Andrew for info that its now standard practice to list outputs before inputs, I had no idea about that! To be fair, I made my new kernel by naively imitating a similar existing kernel (pc2_checks_kernel_mod.F90), and that one had the inputs first. What was the reason for changing to outputs first? (in and then out feels more intuitive to me?)

Cheers!
Mike

[Andrew Coughtrie (andrewcoughtrie)]

The write fields coming first comes mainly from the fact that the nlayers argument is defined by the first field to be passed in and it is usually the writing fields that are being updated and thus need to have looping over the entire column. This is also an historical attempt at consistency across the whole code base where the out/inout arguments should also come first as stated (influenced by that kernel argument order I imagine - though these decisions are long in the past now). The fact that this isn't the case everywhere makes the code harder to read unfortunately. Perhaps it isn't necessary in most subroutines but it is currently safer to do so with kernels, especially if you ever work with 2D fields which have been known to cause headaches when accidently ending up with nlayers=1 due to passing the 2D field first.

This of course isn't an iron clad never break it rule as there are times when what comes first needs to be different (see when-can-i-break-the-rules).

[Steve Mullerworth (stevemullerworth)]

"a) This is exactly the sort of bug that you can catch by testing bit-comparison on different processor decompositions. Why was this key safety-net quietly dropped in the move from um to lfric? We urgently need to implement a bit-comp test for the physics in lfric_apps rose-stem!"

On this point, the option was not "quietly dropped". There was a lot of discussion, particularly with climate scientist strategic heads at the time because historically this ability has been more important for long running climate projections than for NWP. We still need to do work to ensure that we know how to regenerate an existing climate model generated dataset if the requirement arises. This will require remembering what resource was used to run it.

As mentioned before on this discussion group, I'm sure it is feasible to implement bit comparison tests if a suitable physics configuration is required. I believe that, with work to create a reproducible global sum, we could retain bit comparison for different configurations by imposing compute order on the unstructured mesh. This would be for testing only as we would expect there to be an impact on compute performance. It's not something we could find time for now, though.

Just to add to Andy's comment on b), it's not the order of the arguments that have to be right, it's the declarations that have to be right. i.e. real :: my_theta_field(undf_wtheta), where the dummy argument undf_wtheta has to be in the correct place in the argument list. But extracting that information from the kernel code is quite challenging as the declarations can be written in different ways by different people (and there might still be errors in the declaration, even if the argument list is consistently ordered). When development of PSyclone was started, it was decided that the key thing was to base generation on just metadata (the gubbins in the type, public, extends(kernel_type) ::... block) and offer a stub generator for use by people to give them a starting point for writing their kernels.

[MichaelWhitall]

Hi thanks all :) The stub generator thing sounds really handy; I'll try and remember to use that next time I make a new kernel...

But I fear the greatest danger of making errors is not when making a new kernel (when the stub generator can be used to guarantee its all consistent), but when altering the argument list to an existing kernel. It'd be grand if the existing functionality in the stub generator (to autogenerate the argument list and argument declarations consistent with the meta-data) could be deployed in the regular compile? Another reason this would be good is that as it stands, the developer has to specify the list of arguments to a kernel 5 times:

• (a) In the kernel meta-data
• (b) In the kernel's `_code' subroutine argument list
• (c) In the kernel's `_code' subroutine argument declarations
• (d) In the calling _alg code's invocation of the kernel
• (e) In the kernel's doxygen comments describing the argument list

This feels like way too much redundancy / duplication, and too much effort to keep all 5 of these consistent!

Could we add one extra field to the arg_type derived-type structure, so that the kernel meta-data can contain a character string specifying the variable name each argument is given within the kernel? So then instead of writing:

type, public, extends(kernel_type) :: lsp_precfrac_checks_kernel_type
  private
  type(arg_type) :: meta_args(8) = (/                                          &
       arg_type(GH_FIELD, GH_REAL, GH_READ,  W3),     & ! exner_w3
       arg_type(GH_FIELD, GH_REAL, GH_READ,  WTHETA), & ! mv_wth
       arg_type(GH_FIELD, GH_REAL, GH_READ,  WTHETA), & ! ml_wth
       arg_type(GH_FIELD, GH_REAL, GH_READ,  WTHETA), & ! mi_wth
       arg_type(GH_FIELD, GH_REAL, GH_READ,  WTHETA), & ! ms_wth
       arg_type(GH_FIELD, GH_REAL, GH_READ,  WTHETA), & ! mr_wth
       arg_type(GH_FIELD, GH_REAL, GH_READ,  WTHETA), & ! mg_wth
       arg_type(GH_FIELD, GH_REAL, GH_WRITE, WTHETA)  & ! precfrac_wth
       /)
   integer :: operates_on = CELL_COLUMN
contains
  procedure, nopass :: lsp_precfrac_checks_code
end type

(where we've redundantly had to add a comment to remind ourselves which local variable-name each argument corresponds to), we would write:

type, public, extends(kernel_type) :: lsp_precfrac_checks_kernel_type
  private
  type(arg_type) :: meta_args(8) = (/                                          &
       arg_type("exner_w3", GH_FIELD, GH_REAL, GH_READ,  W3),   &
       arg_type("mv_wth", GH_FIELD, GH_REAL, GH_READ,  WTHETA), &
       arg_type("ml_wth", GH_FIELD, GH_REAL, GH_READ,  WTHETA), &
       arg_type("mi_wth", GH_FIELD, GH_REAL, GH_READ,  WTHETA), &
       arg_type("ms_wth", GH_FIELD, GH_REAL, GH_READ,  WTHETA), &
       arg_type("mr_wth", GH_FIELD, GH_REAL, GH_READ,  WTHETA), &
       arg_type("mg_wth", GH_FIELD, GH_REAL, GH_READ,  WTHETA), &
       arg_type("precfrac_wth", GH_FIELD, GH_REAL, GH_WRITE, WTHETA)  &
       /)
   integer :: operates_on = CELL_COLUMN
contains
  procedure, nopass :: lsp_precfrac_checks_code
end type

And then with this information, psyclone would have all the meta-data it needs to autogenerate the kernel _code subroutine's full argument list and argument declarations consistently. That would remove the need for the developer to write (b) and (c) in my above list themselves, reducing the needless duplication and potential for error. Of course, it would still be on the developer to ensure the calling code lists the arguments in the same order as they are listed in the subroutine's argument list, but that's always been the case ;) What do you reckon?

Cheers!
Mike

[Steve Mullerworth (stevemullerworth)]

Compared to normal Fortran, the only additional piece of information to add is the metadata, and the metadata allows checks and code generation that is a benefit as compared with the cost of adding it.

The existing metadata schema was created to support GungHo dynamics. In GungHo, kernels are more compact, and the metadata help to document the kernel. While the schema also works for physics kernels, this was not anticipated. I had envisaged that we would need to do more work to plumb physics into GungHo, but the GHASP project found a way of shoe-horning it in without too much extra work (basically, by reordering the physics field data from UM-style layer-by-layer to GungHo-style column-by-column). This was great for science progress but it has meant that we have a schema for the physics kernels that is not ideal.

There is ongoing work to change the metadata to make it more helpful for the physics kernels, though it won't do precisely what you are asking. The planned change should help with physics kernels that have several dozen fields and a mix of field types. Implementing changes that get the benefit of this change will be a lot of work, so it is hard to know when it will be possible to do it (I don't work closely on this aspect so am not fully aware of plans).