Field vector and array of fields and psy_data

[hiker]

@christophermaynard raised the issue about handling of array of fields. This issue is just to keep track of notes while extending the example. First of all, psydata handles array of fields: it passes the array of fields to the psydata layer, which in turn writes all components of the array. E.g.:

     type(arg_type), dimension(4) :: meta_args =       &
          (/ arg_type(gh_field, gh_real, gh_inc,  w0), &
             arg_type(gh_field, gh_real, gh_read, w0), &
             arg_type(gh_field*3, gh_real, gh_read, w0), &
             arg_type(gh_scalar, gh_logical, gh_read)  &
           /)
...
      TYPE(field_type), intent(in) :: field1, field2, chi(3)
...
      chi_proxy(1) = chi(1)%get_proxy()
      chi_proxy(2) = chi(2)%get_proxy()
      chi_proxy(3) = chi(3)%get_proxy()
...
      CALL extract_psy_data%PreDeclareVariable("chi", chi)

And it's then up to the psydata layer to handle this. For extraction (and, without testing, all other psydata transformations) it will loop across all fields in the array and calls the corresponding function for a field. For example, in extraction:

        type(field_type), dimension(:), intent(in)        :: value
...
        ! Provide each dimension of the vector as an individual 1d array.
        do i = 1, size(value, 1)
            value_proxy = value(i)%get_proxy()
            write(number, '("%",i0)') i
            call this%ProvideVariable(name//trim(number), value_proxy%data)
        enddo

So a 3-element array of fields chi, will add chi%1, chi%2, and chi%3 to the output file ('%' is indeed a valid component of a name in NetCDF). In the output file:

dimensions:
        chi%1dim%1 = 539 ;
        chi%2dim%1 = 539 ;
        chi%3dim%1 = 539 ;
...
variables:
        double chi%1(chi%1dim%1) ;
        double chi%2(chi%2dim%1) ;
        double chi%3(chi%3dim%1) ;
 chi%1 = 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
...
 chi%2 = 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
...
 chi%3 = 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 

Standalone will have the same information (just no names to distinguish them, it is all order dependent).

I believe once we get the kernel extraction working, this might need some improvements, e.g. it might be more convenient to also store the number of fields in an array - atm the number is determined at runtime and would therefore be required, but it coulso be that the driver can get the right number from the kernel meta-data.

[hiker]

While there is a field_vector_type in LFRic, this type appears not to be used (or supported??) in PSyclone. The typical usage pattern is:

    select type (x)
    type is (field_vector_type)
      select type (y)
      type is (field_vector_type)
        ! Get pointers to x and y abstract types
        x_vec => x%get_field(1)
        y_vec => y%get_field(1)

        ! Compute M3*p + H[ D{\theta*}*p ]
        if ( element_order == 0 ) then
          Helmholtz_operator => get_helmholtz_operator(self%level)
          call invoke( name="apply_H",&
                       apply_helmholtz_operator_kernel_type(y_vec,   &
                                                            x_vec,   &
                                                            one_int, &
                                                            Helmholtz_operator, &
                                                            lam_mesh) )

@christophermaynard, @TeranIvy, @arporter, is it correct to state: The individual vectors of a field_vector_type are extracted and then individually provided as fields. The gh_field*X notation appears to be restricted to array of fields only, and there is no direct support of field_vector_type.

The examples/tests in PSyclone seems to use all the same approach, e.g.:

  type :: some_type
     type(r_solver_field_vector_type) :: vec_type(10)
   contains
     procedure, public :: my_sub
  end type some_type

  contains

  subroutine my_sub(self, x, m1, m2)
    class(some_type), intent(inout) :: self
...
    call invoke(testkern_type(a, x%vector(1), self%vec_type(1)%vector(1), m1, m2))

If this is correct, we might need to work on our documentation a bit (dynamo.rst), which frequently mentions 'vector', while it very likely means just an array of fields. Same applies to my psydata code - I think I called functions even DeclareFieldVector etc, while de-facto these should be some kind of FieldArrays.

[arporter]

I agree - we do need to update the documentation because we've always called e.g. gh_field*3 a vector. This does beg the question what the vector type is and, as you say, whether we'll need to support it.

[hiker]

These two concepts seems indeed to be very similar (main difference as far as I can tell: the FieldVector provides linear algebra functions to be applied to all its individual vectors. If it needs to be supported, it should be very very similar to what is already there.

From my point of view: PSyData will be fine with it (with some minor changes), and it would be trivial to add support if PSyclone ever explicitly supports FieldVector.

[rupertford]

The problem with field_vector_type is that it is abstract and therefore it is not possible for PSyclone to determine the datatype of anything declared as a field_vector_type in the algorithm layer. The workaround is for the scientists to specify pointers with the concrete type which they associate within the appropriate select (see x_vec and y_vec in your example). PSyclone is then able to determine the datatype of the arguments. So I'm not sure it is right to say that PSyclone does not support field_vector_type, in the sense that there is algorithm code that uses it, it just does not make use of it in the PSy layer.

[hiker]

The problem with field_vector_type is that it is abstract and therefore it is not possible for PSyclone to determine the datatype of anything declared as a field_vector_type in the algorithm layer.

Technically field_vector_type is not abstract, but mostly LFRic is indeed using the abstract_vector_type. And it took me a while to see where the problem with the abstract type is - I assume it is to be able to properly declare the proxy type?

So does that mean that if I have a field_vector_type say a (not abstract), I can pass this into an a kernel that expects a gh_Field*3 (or whatever)? And PSyclone will create say a1_proxy = a%vector(1)%get_proxy(), a2_proxy=a%vector(2)%get_proxy() a3_proxy=..., and then pass the a_proxy*%data all as fields into the kernel call? Guess I can just try that .... and I get:

Generation Error: The metadata for argument 'chi_fv' in kernel 'testkern_w0_code' specifies that this is a real field, however it is declared as a 'field_vector_type' in the algorithm code.

The workaround is for the scientists to specify pointers with the concrete type which they associate within the appropriate select (see x_vec and y_vec in your example). PSyclone is then able to determine the datatype of the arguments. So I'm not sure it is right to say that PSyclone does not support field_vector_type, in the sense that there is algorithm code that uses it, it just does not make use of it in the PSy layer.

Thanks for the clarification, this makes sense then.

[rupertford]

You're quite right @hiker, I was think of the abstract_vector_type. The field_vector_type is used by the scientist in the algorithm layer and they must dereference in the algorithm layer to individual fields to use in PSyclone.