Search code examples
matrixjuliabinaryfiles

Reading unformatted FORTRAN files into multidimensional arrays using Julia

I'm just getting started with Julia, and I'm trying to read an unformatted FORTRAN file and store the data in arrays that are shaped in a particular way. I'm not sure how to accomplish this using Julia.

I've found the Julia package FortranFiles, which provides a direct way to read unformatted FORTRAN files using Julia. The file I'm trying to read looks like:

1 integer:
[nzones]

nzones*3 integers (brackets indicate one record):
[idim1,jdim1,kdim1,idim2,jdim2,kdim2,...,
idim_nzones,jdim_nzones,kdim_nzones]

series of nzones datasets:
[xvalues1,yvalues1,zvalues1](floating point values) for 1st zone
[xvalues1,yvalues1,zvalues1](floating point values) for 2nd zone
...,
[xvalues1,yvalues1,zvalues1](floating point values) for last zone

where the first line represents the number of zones and the lines that follow represent a grid dimension in each i, j, and k directions. Following these first two records are the x, y, and z coordinates, which are Float64s, for each i, j, and k point in a zone, and I would like to shape the arrays as x(1:im,1:jm,1:km,m), y(1:im,1:jm,1:km,m), and z(1:im,1:jm,1:km,m) where im, jm, and km are the imax,jmax, and kmax extents listed for each zone. Here's what I have so far:

using FortranFiles

fname = "my_file"
fid = FortranFile(fname)

@fread fid nblks::Int32
@fread fid ni::(Int32,nblks) nj::(Int32,nblks) nk::(Int32,nblks)

Here's where I'm getting hung up. For each zone I have x, y, and z coordinate arrays which should all be rank 4 arrays. For the x array, I want to store all of the x coordinates where x[1,1,1,1] refers to an x coordinate value at i=1, j=1, k=1, and zone =1 and x[end, end, end, end] refers to an x coordinate value at i = imax, j=jmax, k=kmax, and for the last zone listed (i.,e. zone = nblks). Then I want to create similar arrays for the y and z coordinate values.

Something like:

for m = 1:nblks
    im = ni[m]
    jm = nj[m]
    km = nk[m]

    @fread fid x::(Float64,im,jm,km,m) y::(Float64,im,jm,km,m) z::(Float64,im,jm,km,m)
end

However, I get a FortranFilesError: attempting to read beyond record end when trying this approach.

Solution

  • It appears that my issue is somewhat related to how Julia reads unformatted binary data, which is different from how FORTRAN's read works on the same data.

    In FORTRAN, I could do something like:

    integer, dimension (:), allocatable :: idim, jdim, kdim
integer :: nblks, fid, ios

fid = 10
open(unit=fid,form='unformatted', file='my_file',status='old',iostat=ios)

if( ios /= 0 ) then
    write(*,*) '*** Error reading file ***'
    stop
end if

read(fid) nblks

allocate( idim(nblks), jdim(nblks), kdim(nblks) )

read(fid) ( idim(m), jdim(m), kdim(m), m = 1, nblks )

close(fid)
...

    However in Julia, I need to keep track of the file pointer's position, and realize that each record is preceded and followed by a 4-byte integer. I haven't been able to find a way to read each zone's i, j, & k extents directly into three separate arrays like can be done in FORTRAN (since the record is probably parsed line by line), but an alternative in Julia is to just read the entire record into a single nblk*3 element vector, and then reshape this vector afterwards:

    fid = open("my_file")

skip(fid,4)

nblks = read(fid,Int32)

skip(fid,8)

dims = Array{Int32}(undef,3*nblks)

read!(fid,dims)

ni, nj, nk = [Array{Int32}(undef,nblks) for i in 1:3]

for m in 1:nblks
    ni[m] = dims[3*m-2]
    nj[m] = dims[3*m-1]
    nk[m] = dims[3*m]
end