Global subroutine to import fields from fld files (gfldr)

[vbaconnet]

As the title says. Introduces a global subroutine to handle reading + interpolation of fld files, similar in essence to the Nek5000 gfldr, but on steroids.

Fixes #2284 and #2268

The global subroutine is named import_fields. Below are some examples of how the function can be used.

! ---- Say that u,v,w,p,t, etc are pointers to whatever fields
type(field_t), pointer :: u,v,w,p,t,s1,s2,my_custom_field
type(field_list_t) :: my_field_list
! ----

! Case 1. Import velocity and pressure (no interpolation, only potentially space-to-space interpolation)
call import_fields("field0.f00001", u=u, v=v, w=w, p=p)

! Case 2. Import velocity and pressure (with interpolation)
call import_fields("field0.f00001", u=u, v=v, w=w, p=p, interpolation=.true., tolerance=1e-9)

! Case 2 (bis). Import velocity and pressure (with interpolation and specified mesh file)
call import_fields("field0.f00001", u=u, v=v, w=w, p=p, interpolation=.true., tolerance=1e-9, mesh_fname="field0.f00000")

! Case 3. Import the two first velocity arrays, but take the w velocity from the temperature array.
call import_fields("field0.f00001", u=u, v=v, &
     t=w) ! this means "the temperature array in "field0.f00001" should be loaded into the field w" 

! Case 4. "field0.f00001" has velocity fields which I want to load into random fields.
call import_fields("field0.f00001", u=s1, v=s2, w=my_custom_field) 

! Case 5. "field0.f00001" contains temperature, s01, s02. 
call my_field_list%init(2)
call my_field_list%assign(1, s1)
call my_field_list%assign(2, s2)
call import_fields("field0.f00001", t=my_custom_field, s_target_list=my_field_list) 
call my_field_list%free()

! Case 5 (bis). Same as case 5 but load temperature from the scalar list instead (equivalent)
call my_field_list%init(3)
call my_field_list%assign(1, my_custom_field)
call my_field_list%assign(2, s1)
call my_field_list%assign(3, s2)
call import_fields("field0.f00001", s_target_list=my_field_list, s_index_list=(/0,1,2/)) ! index 0 loads the temperature field from the fld file.
call my_field_list%free()

! Case 6. "field0.f00001" has arrays s01, s02, s03 which I want to load into u,v,w. 
call my_field_list%init(3)
call my_field_list%assign(1, u)
call my_field_list%assign(2, v)
call my_field_list%assign(3, w)
call import_fields("field0.f00001", s_target_list=my_field_list) 
call my_field_list%free()

• As you can see this routine makes use of optional arguments to precisely select which fields should be loaded, and from which array in the fld file. I guess we haven't used things like this in the past so let me know what you think. This is very much inspired by the pynekread/pynekwrite functions from pySEMTools!
• I think this is quite powerful and can be used in a lot of places in the code. I have already changed the initial conditions and the sponge source term accordingly.
• The name import_fields is very much open for discussion, but I would suggest to keep it short, as you can see the argument list can grow quickly. Same goes for the dummy arguments s_target_list and s_index_list.
• Fixes Field IC doesn't work with a single scalar if it is not temperature #2284 : introduces a new parameter target_index (maybe should be source_index) to the scalar initial condition, that allows the user to bypass the naming conventions and select which scalar/temperature array they want to import their data from. target_index = 0 load the temperature array in the fld file. So you effectively don't need your field to be named temperature in the registry to load the data from the temperature array.

More details

• import_fields is implemented in its own file because it was a mess with all different imports: file, fld_file_data, global_interpolation, etc. But I guess that's a choice that can be discussed
• import_fields is a basically fancy wrapper ("one-call wonder") that takes care of all the machinery of file names, input arguments, data movement (memcpy) etc. The interpolation operations are offloaded to the new fld_file_data_import_fields subroutine.

About fld_file_data_import_fields

• This subroutine does mesh-to-mesh or space-to-space interpolation, but it's essentially a glorified copy of fld arrays to fields.
• It does NOT do any data movement, that is handled by the global subroutine import_fields.

One thing that was bugging me is if, and where data movement should be done. When an fld file is read and all the data is loaded into the fld_file_data object, there is no device memcpy. Since the interpolation is done on GPUs, a memcpy would be required after an fld_file read:

call f%read(fld_data)
! ---> Need a memcpy of x,y,z and u here! <---
call fld_data%import_fields(u=u, interpolate=.true.)

So, should fld_file_data_import_fields take care of data movement for you? Well thankfully we have a note about data movement in our documentation:

It is the programmers' responsibility to make sure that device arrays are kept in sync with the associated host array. Neko does not perform any implicit data movement.

Therefore this routine is kept lightweight and focuses only on interpolation capabilities. Instead, the data movement aspect is taken care of by the global import_fields subroutine. Because it takes the file name as input, there is no access to the fld_file_data object and no possibility to handle the memcpy manually anyways.

Compiling

• GNU fortran
• CUDA
• HIP/Cray

Testing

• GNU fortan
• CUDA
• HIP/Cray

[njansson]

Nice!

[timofeymukha]

Awesome stuff!! Only glanced through the code. But dont call it global, this is a totally normal subroutine living in a module. The word global scares people :-).

[njansson]

For the future, would be very nice if one can also support hdf5 files ;)

[timfelle]

For the future, would be very nice if one can also support hdf5 files ;)

I have been working a bit with HDF5 now, and it should be quite simple to extend, atleast to be functional. But we would need to decide on some standardized path naming, Which should probably be informed by however a vtkhdf5 neko export would look like.

[adperezm]

Regarding the name, I would just call it read_fields or so :)

[vbaconnet]

@adperezm replying to your PR #2327 about vector lists here:

I think it would be nice to add names to the vectors. [...] I would like that functionality to accept a list of keys to read, and since fld_file_data_t uses vector_t, it would be nice if one could search a vector by name. This would also help to have a unique interface with hdf5, I believe.

I do not like that these searches need to always traverse the list. But I suppose that having hashtables with these types is overkill, considering that we never have crazy many entries.

Yea I think this would be definitely great for an hdf5 interface.

For the fld files my opinion is that we shouldn't spend more time on extending fld_file_data_t with hash tables, if we're already trying to move away from that format? All of those vectors are hardcoded to mean something specific anyways, and I don't know if passing a list of scalar indices (/0,4,5/) is fundamentally different from (/"temperature", "s04", "s05"/), except being more verbose? Unless we foresee to actually add metadata to the nek5000 file with field labels, which somehow would end up as vector names after being read, but that seems like a mess. But maybe I don't see the angle here... :)

And sorry I don't understand what you mean about the searches traversing the lists, for import_fields there is no search since the user specifies the exact fields they need, but I probably missed your point :D

[adperezm]

@adperezm replying to your PR #2327 about vector lists here:

I think it would be nice to add names to the vectors. [...] I would like that functionality to accept a list of keys to read, and since fld_file_data_t uses vector_t, it would be nice if one could search a vector by name. This would also help to have a unique interface with hdf5, I believe.
I do not like that these searches need to always traverse the list. But I suppose that having hashtables with these types is overkill, considering that we never have crazy many entries.

Yea I think this would be definitely great for an hdf5 interface.

For the fld files my opinion is that we shouldn't spend more time on extending fld_file_data_t with hash tables, if we're already trying to move away from that format? All of those vectors are hardcoded to mean something specific anyways, and I don't know if passing a list of scalar indices (/0,4,5/) is fundamentally different from (/"temperature", "s04", "s05"/), except being more verbose? Unless we foresee to actually add metadata to the nek5000 file with field labels, which somehow would end up as vector names after being read, but that seems like a mess. But maybe I don't see the angle here... :)

And sorry I don't understand what you mean about the searches traversing the lists, for import_fields there is no search since the user specifies the exact fields they need, but I probably missed your point :D

No problem :), The part about the search was not specific to this PR. Just in general, when we search for something by name in a list, we traverse, so there is a bunch of implicit nested loops there. But this does not affect this PR at all.

The reason I was thinking about extending this function you provided here, was because in my mind, something like this would be nice (highjacking your example):

! Case 5. "field0.f00001" contains temperature, s01, s02.

! Define a list with sources
character(:), allocatable :: names(:)
names = [ character(:) :: "t", "s1", "s2" ]

! Define a list with outputs (potentially also replace with field_arrays_t that also initialize the data)
call my_field_list%init(3)
call my_field_list%assign(1, my_custom_field)
call my_field_list%assign(2, s1)
call my_field_list%assign(3, s2)

! Do the import
call import_fields("field0.f00001", output=my_field_list, from=names)

And that would be doable if one has names for the vectors, even if it is hard-coded. I think it would nicely change things inside the current code from:

 ! Evaluate all the fields
       if (present(u)) call global_interp%evaluate(u%x(:,1,1,1), this%u%x, &
               on_host=.false.)
       if (present(v)) call global_interp%evaluate(v%x(:,1,1,1), this%v%x, &
               on_host=.false.)
       if (present(w)) call global_interp%evaluate(w%x(:,1,1,1), this%w%x, &
               on_host=.false.)
       if (present(p)) call global_interp%evaluate(p%x(:,1,1,1), this%p%x, &
               on_host=.false.)
       if (present(t)) call global_interp%evaluate(t%x(:,1,1,1), this%t%x, &
               on_host=.false.) 
       if (present(s_target_list)) then

          ! If the index list exists, use it as a "mask"
          if (present(s_index_list)) then
             do i = 1, size(s_index_list)
                ! Take care that if we set i=0 we want temperature
                if (s_index_list(i) .eq. 0) then
                   call global_interp%evaluate(s_target_list%x(i), &
                           this%t%x, on_host=.false.)
                else
                   call global_interp%evaluate(s_target_list%x(i), &
                           this%s(s_index_list(i))%x, on_host=.false.)
                end if
             end do

To something like (not really proper fortran):

 ! Evaluate all the fields
do i = 1, n
 call global_interp%evaluate(output_field_list%get_by_index(i)%x(:,1,1,1), &
                                                   this%vector_list%get_by_name(names(i))%x, &
                                                  on_host=.false.)
end do

If one then uses this import_fields routine but gives a hdf5 as file, then one could search inside for keys as well.

This is basically what I do in pysemtools with the function read_data (https://github.com/ExtremeFLOW/pySEMTools/blob/59350bd439c79f8765fb7560388447bbda63cd2d/pysemtools/io/wrappers.py#L70)

I have found very useful to use the same simple interface for any type of file I want to read.

But in the end, I did not say we should change this :). It does come to preferences. I would like to always output a list of fields and input a list of keys, then I do not need to play around with optional arguments in the subroutine, but this does not mean everyone else likes that.

[vbaconnet]

To something like (not really proper fortran):

 ! Evaluate all the fields
do i = 1, n
 call global_interp%evaluate(output_field_list%get_by_index(i)%x(:,1,1,1), &
                                                   this%vector_list%get_by_name(names(i))%x, &
                                                  on_host=.false.)
end do

oooh I see now. Didn't think that far, sorry.

I have found very useful to use the same simple interface for any type of file I want to read.

Yea definitely agree, once this moves further to accept hdf5 it will for sure be needed

[timofeymukha]

Another alternative could be load_fields. I like it more than read_ since we potentially do extra stuff. But import_ is also not that bad!

[vbaconnet]

I do like import_fields as well.. but let's finalize this in the dev meeting tomorrow so this can move forward to merging!

[vbaconnet]

Unless there are strong opinions I then suggest we jsut go ahead with that name. I wrote that the name was up for discussion in case it was divisive but it doesn't look that way :)

[timofeymukha]

All good, just some tiny stuff.