user3394147
user3394147

Reputation: 107

Accessing data stored in complex arrays

I'm writing some fortran code where I'm having difficulty visualizing how the data in a hierarchy of arrays is stored, which is causing me troubles with how to manipulate a subset of this data. I've read-in some data stored in an unformatted binary PLOT3D file. The data format looks like the following:

  1. 1 int (nblocks): number of blocks (or zones) in a computational mesh
  2. 4 x nblocks ints: (ni(m), nj(m), nk(m), nvars(m)): number of i, j, and k points in each block along with a number of solution variables in each block.
  3. loop over blocks (m)
  4. ni(m) x nj(m) x nk(m)x nvars(m) floats (q(ni,nj,nk,nv,m)): solution variables at each i, j, and k points for each block.
  5. end loop over blocks
program my_program

use iso_fortran_env

implicit none

character(*), parameter :: soln_file = 'my_file_name.q'

integer :: nblks, io_stat, imax, jmax, kmax, nv, m

integer, dimension (:), allocatable :: ni, nj, nk, nvars

real(real64), dimension (:,:,:,:), allocatable :: qq
real(real64), dimension (:,:,:,:,:), allocatable :: q

open(unit=10, form='unformatted', file=soln_file, status='old', iostat=io_stat)

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

read(10) nblks

allocate( ni(nblks), nj(nblks), nk(nblks) )

read(10) ( ni(m), nj(m), nk(m), nvars(m), m = 1, nblks )

imax = maxval(ni)
jmax = maxval(nj)
kmax = maxval(nk)
nv = maxval(nvars)

allocate( q(imax,jmax,kmax,nv,nblks) )

do m = 1, nblks
    allocate( qq(ni(m),nj(m),nk(m),nvars(m)) )

    read(10) qq(ni(m),nj(m),nk(m),nvars(m))

    q(1:ni(m),1:nj(m),1:nk(m),1:nvars(m),m) = qq

    deallocate( qq )
end do

close(10)

deallocate( ni, nj, nk, nvars, q )

stop

end program my_program

In my case, I'm interested in extracting a subset, or just modifying the values in place, of a single solution variable at all of the points in every block. It seems like I should end up with the solution values for that variable at each i, j, and k point for every block in a 3 x nblocks array.

However, I get a seg fault at run time which indicates to me that I'm not sizing the array correctly.

Upvotes: 2

Views: 110

Answers (1)

PTRK
PTRK

Reputation: 929

In your block of code

do m = 1, nblks
    allocate( qq(ni(m),nj(m),nk(m),nvars(m)) )

    read(10) qq(ni(m),nj(m),nk(m),nvars(m))

    q(1:ni(m),1:nj(m),1:nk(m),1:nvars(m),m) = qq

    deallocate( qq )
end do

You're only reading 4 values pour qq, instead of ni(m) x nj(m) x nk(m) x nvars(m). Then you're trying to copy theses values into q wich is not conformant.

It seems to me that the correct loop should be

do m = 1, nblks
    allocate( qq(ni(m),nj(m),nk(m),nvars(m)) )

    do iv = 1 ,nvars(m)
      do ik = 1 ,nk(m)
        do ij = 1 ,nj(m)
          do ii = 1 ,ni(m)
             read(10) qq(ii,ij,ik,iv)
          end do
        end do
      end do
    end do

    q(1:ni(m),1:nj(m),1:nk(m),1:nvars(m),m) = qq

    deallocate( qq )
end do

I'm not sure this will solve your problem, but there's something fishy here.

Upvotes: 1

Related Questions