readgrav_gr.py

# This script read data from grav_ files and save to .npy
# Fortran unformatted files seem to be composed with several blocks of different datatype.
# Each block looks like this:
# |(int) Size of DATA || (dtype) DATA ||(int) Size of DATA |

import numpy as np
from struct import *
import time
import sys

#t1=time.process_time()
#t1=time.process_time()
#simulation details
levelmin=int(sys.argv[3])
levelmax=levelmin
twotondim=8
G=int(sys.argv[4])

#files to be read
path=sys.argv[1]
#outpath=sys.argv[2]
fileroot='grav_'+sys.argv[1][-6:-1]+'.out'
numfiles=int(sys.argv[2])
firstcount=1
outfile=fileroot

#internal variable
datlis=[phi,f1,f2,f3,gr5,gr6,gr7,gr8,gr9,gr10,rho,gr_m1,gr_m2,gr_m3,gr_m4]=[[],[],[],[],[],[],[],[],[],[],[],[],[],[],[]]
print(len(datlis))
Ngrid=0
for i in range(firstcount, firstcount+numfiles):
    infile=fileroot+str(i).zfill(5)
    with open(path+infile,'rb') as f:
        content = f.read()
#    print('Reading '+path+infile+'  : '+str(time.process_time()-t1)+' s')
    pmin=0
    pmax=48

    ############ 'header info'############
    info=unpack('i'*3*4,content[pmin:pmax])
    [ncpu,ndim,nlevelmax,nboundary]=[info[1],info[4],info[7],info[10]]

    ############ levels ############
    for ilevel in range(levelmin,levelmax+1):
        ############ nboundary+ncpu ############
        for ibound in range(1,nboundary+ncpu+1):
            pmin0=pmax
            pmax0=pmin0+4*3*2
            info=unpack('i'*3*2,content[pmin0:pmax0])
            [currlevel,ncache]=[info[1],info[4]]
            Ngrid+=ncache*twotondim
            #print pmin0, pmax0, currlevel, ncache
            if ncache==0:
                pmax=pmax0
                continue
            ############ twotondim ############
            for ind in range(1,twotondim+1):
                j=0
                ############ parameters ############
                for N in ([1,2,3,4,5,6,7,8,9,10,11,12,13,14,15]):
                    ipar=datlis[j]
                    pmin=pmax0+(8*N-4)+(N-1)*8*ncache
                    pmax=pmin+ncache*8
                    #print(pmin,pmax,ncache)
                    info=unpack('d'*ncache,content[pmin:pmax])
                    for floatelem in info:
                        ipar.append(floatelem)
                    j=j+1
                    #pmax=pmax+4
                    #print ind, pmin, pmax, ncache
                pmax0=pmax+4
            pmax=pmax0
    f.close()

#print('Transposing: '+str(time.process_time()-t1)+' s')
datlis=np.transpose(datlis)
print(len(datlis))
#print('Deleting repeating elements: '+str(time.process_time()-t1)+' s')
#datlis=np.vstack({tuple(row) for row in datlis})
#print('Saving to '+path+outfile+'  : '+str(time.process_time()-t1)+' s')
#np.save(path+outfile,datlis)

#print('Mapping to tuple: '+str(time.process_time()-t1)+' s')

datlis=list(map(tuple,datlis))
print(len(datlis))

#print('Finding set of values: '+str(time.process_time()-t1)+' s')
datlis=set(datlis)
print(len(datlis))

#print('Mapping back to array: '+str(time.process_time()-t1)+' s')
datlis=list(map(tuple,datlis))
print(len(datlis))

data=np.zeros((G**3,15))
for i in range(len(data)):
    print(i,len(data),len(datlis))
    data[i][0]=datlis[i][0]
    data[i][1]=datlis[i][1]
    data[i][2]=datlis[i][2]
    data[i][3]=datlis[i][3]
    data[i][4]=datlis[i][4]
    data[i][5]=datlis[i][5]
    data[i][6]=datlis[i][6]
    data[i][7]=datlis[i][7]
    data[i][8]=datlis[i][8]
    data[i][9]=datlis[i][9]
    data[i][10]=datlis[i][10]
    data[i][11]=datlis[i][11]
    data[i][12]=datlis[i][12]
    data[i][13]=datlis[i][13]
    data[i][14]=datlis[i][14]

#print('Saving to '+path+outfile+'  : '+str(time.process_time()-t1)+' s')
np.save(path+outfile,data)

#print('Done : '+str(time.process_time()-t1)+' s, N = %d'%(len(datlis)))
sys.exit()