228 lines
8.4 KiB
Python
228 lines
8.4 KiB
Python
# functinos for post-processing output data in h5 file
|
|
|
|
import numpy
|
|
import h5py
|
|
|
|
|
|
def get_data_from_h5(fpath, fpath_grid, dataset_name, nr_glob, nt_glob, np_glob, ndiv_r, ndiv_t, ndiv_p, verbose=False):
|
|
"""
|
|
Output arrays has one overlapping layer between adjusted subdomains.
|
|
This function will skip those layers and reconstruct the entire grid dataset.
|
|
|
|
fpath: path to field data file
|
|
fpaht_grid: path to grid data file
|
|
dataset_name: name of dataset in h5 file
|
|
nr_glob: number of grid points in r direction
|
|
nt_glob: number of grid points in t direction
|
|
np_glob: number of grid points in p direction
|
|
ndiv_r: number of subdomains in r direction
|
|
ndiv_t: number of subdomains in t direction
|
|
ndiv_p: number of subdomains in p direction
|
|
verbose: print out information
|
|
|
|
"""
|
|
# total number of subdomains
|
|
n_sub = ndiv_r*ndiv_t*ndiv_p
|
|
|
|
# total points on each direction with overlap
|
|
nr_total_glob = nr_glob + ndiv_r - 1
|
|
nt_total_glob = nt_glob + ndiv_t - 1
|
|
np_total_glob = np_glob + ndiv_p - 1
|
|
|
|
# prepare a 3D array to store the data
|
|
data_glob = numpy.zeros((nr_glob,nt_glob,np_glob), dtype=numpy.float64)
|
|
grid_glob_r = numpy.zeros((nr_glob,nt_glob,np_glob), dtype=numpy.float64)
|
|
grid_glob_t = numpy.zeros((nr_glob,nt_glob,np_glob), dtype=numpy.float64)
|
|
grid_glob_p = numpy.zeros((nr_glob,nt_glob,np_glob), dtype=numpy.float64)
|
|
|
|
# open grid data file
|
|
fgrid = h5py.File(fpath_grid, 'r')
|
|
|
|
# open field data file
|
|
fdata = h5py.File(fpath, 'r')
|
|
|
|
#try:
|
|
# load data data by each subdomain
|
|
|
|
# offset
|
|
offset = 0
|
|
|
|
for ir_sub in range(ndiv_r):
|
|
for it_sub in range(ndiv_t):
|
|
for ip_sub in range(ndiv_p):
|
|
|
|
# number of data point for this sub domain
|
|
nr_sub = nr_glob//ndiv_r
|
|
nt_sub = nt_glob//ndiv_t
|
|
np_sub = np_glob//ndiv_p
|
|
|
|
# offset for each direction
|
|
offset_r = ir_sub*nr_sub
|
|
offset_t = it_sub*nt_sub
|
|
offset_p = ip_sub*np_sub
|
|
|
|
# add modulus to the last subdomains
|
|
if ir_sub == ndiv_r-1:
|
|
nr_sub += nr_glob%ndiv_r
|
|
if it_sub == ndiv_t-1:
|
|
nt_sub += nt_glob%ndiv_t
|
|
if ip_sub == ndiv_p-1:
|
|
np_sub += np_glob%ndiv_p
|
|
|
|
# add overlap layer if this subdomain is not the last one for each direction
|
|
if ir_sub != ndiv_r-1:
|
|
nr_sub += 1
|
|
if it_sub != ndiv_t-1:
|
|
nt_sub += 1
|
|
if ip_sub != ndiv_p-1:
|
|
np_sub += 1
|
|
|
|
# number of data point for this sub domain
|
|
n_points_total_sub = nr_sub*nt_sub*np_sub
|
|
|
|
# load data
|
|
data_sub = fdata[dataset_name][offset:offset+n_points_total_sub]
|
|
data_sub_p = fgrid["/Mesh/node_coords_p"][offset:offset+n_points_total_sub]
|
|
data_sub_t = fgrid["/Mesh/node_coords_t"][offset:offset+n_points_total_sub]
|
|
data_sub_r = fgrid["/Mesh/node_coords_r"][offset:offset+n_points_total_sub]
|
|
|
|
# reshape data
|
|
data_sub = data_sub.reshape(nr_sub, nt_sub, np_sub)
|
|
data_sub_p = data_sub_p.reshape(nr_sub, nt_sub, np_sub)
|
|
data_sub_t = data_sub_t.reshape(nr_sub, nt_sub, np_sub)
|
|
data_sub_r = data_sub_r.reshape(nr_sub, nt_sub, np_sub)
|
|
|
|
|
|
# put those data in global 3d array
|
|
data_glob[offset_r:offset_r+nr_sub, offset_t:offset_t+nt_sub, offset_p:offset_p+np_sub] = data_sub
|
|
grid_glob_p[offset_r:offset_r+nr_sub, offset_t:offset_t+nt_sub, offset_p:offset_p+np_sub] = data_sub_p
|
|
grid_glob_t[offset_r:offset_r+nr_sub, offset_t:offset_t+nt_sub, offset_p:offset_p+np_sub] = data_sub_t
|
|
grid_glob_r[offset_r:offset_r+nr_sub, offset_t:offset_t+nt_sub, offset_p:offset_p+np_sub] = data_sub_r
|
|
|
|
# update offset
|
|
offset += n_points_total_sub
|
|
|
|
|
|
return data_glob, grid_glob_r, grid_glob_t, grid_glob_p
|
|
|
|
#except:
|
|
# fdata.close()
|
|
# fgrid.close()
|
|
# print("error occured while reading file.")
|
|
|
|
# return -1
|
|
|
|
|
|
# output arrays has one overlapping layer between adjusted subdomains
|
|
# this function will skip those layers and reconstruct the entire grid dataset
|
|
def get_data_from_ascii(fpath, fpath_grid, nr_glob, nt_glob, np_glob, ndiv_r, ndiv_t, ndiv_p, verbose=False):
|
|
"""
|
|
|
|
fpath: path to ascii data file
|
|
fpath_grid: path to ascii grid data file
|
|
nr_glob: number of grid points in r direction
|
|
nt_glob: number of grid points in t direction
|
|
np_glob: number of grid points in p direction
|
|
ndiv_r: number of subdomains in r direction
|
|
ndiv_t: number of subdomains in t direction
|
|
ndiv_p: number of subdomains in p direction
|
|
verbose: print out information
|
|
|
|
"""
|
|
# total number of subdomains
|
|
n_sub = ndiv_r*ndiv_t*ndiv_p
|
|
|
|
# total points on each direction with overlap
|
|
nr_total_glob = nr_glob + ndiv_r - 1
|
|
nt_total_glob = nt_glob + ndiv_t - 1
|
|
np_total_glob = np_glob + ndiv_p - 1
|
|
|
|
# prepare a 3D array to store the data
|
|
data_glob = numpy.zeros((nr_glob,nt_glob,np_glob), dtype=numpy.float64)
|
|
grid_glob_r = numpy.zeros((nr_glob,nt_glob,np_glob), dtype=numpy.float64)
|
|
grid_glob_t = numpy.zeros((nr_glob,nt_glob,np_glob), dtype=numpy.float64)
|
|
grid_glob_p = numpy.zeros((nr_glob,nt_glob,np_glob), dtype=numpy.float64)
|
|
|
|
# read data
|
|
data_tmp = numpy.loadtxt(fpath)
|
|
grid_tmp = numpy.loadtxt(fpath_grid)
|
|
|
|
# load data data by each subdomain
|
|
|
|
# offset
|
|
offset = 0
|
|
|
|
for ir_sub in range(ndiv_r):
|
|
for it_sub in range(ndiv_t):
|
|
for ip_sub in range(ndiv_p):
|
|
|
|
# number of data point for this sub domain
|
|
nr_sub = nr_glob//ndiv_r
|
|
nt_sub = nt_glob//ndiv_t
|
|
np_sub = np_glob//ndiv_p
|
|
|
|
# offset for each direction
|
|
offset_r = ir_sub*nr_sub
|
|
offset_t = it_sub*nt_sub
|
|
offset_p = ip_sub*np_sub
|
|
|
|
# add modulus to the last subdomains
|
|
if ir_sub == ndiv_r-1:
|
|
nr_sub += nr_glob%ndiv_r
|
|
if it_sub == ndiv_t-1:
|
|
nt_sub += nt_glob%ndiv_t
|
|
if ip_sub == ndiv_p-1:
|
|
np_sub += np_glob%ndiv_p
|
|
|
|
# add overlap layer if this subdomain is not the last one for each direction
|
|
if ir_sub != ndiv_r-1:
|
|
nr_sub += 1
|
|
if it_sub != ndiv_t-1:
|
|
nt_sub += 1
|
|
if ip_sub != ndiv_p-1:
|
|
np_sub += 1
|
|
|
|
# number of data point for this sub domain
|
|
n_points_total_sub = nr_sub*nt_sub*np_sub
|
|
|
|
# load data
|
|
data_sub = data_tmp[offset:offset+n_points_total_sub]
|
|
grid_sub_p = grid_tmp[offset:offset+n_points_total_sub,0]
|
|
grid_sub_t = grid_tmp[offset:offset+n_points_total_sub,1]
|
|
grid_sub_r = grid_tmp[offset:offset+n_points_total_sub,2]
|
|
|
|
# reshape data
|
|
data_sub = data_sub.reshape(nr_sub, nt_sub, np_sub)
|
|
grid_sub_p = grid_sub_p.reshape(nr_sub, nt_sub, np_sub)
|
|
grid_sub_t = grid_sub_t.reshape(nr_sub, nt_sub, np_sub)
|
|
grid_sub_r = grid_sub_r.reshape(nr_sub, nt_sub, np_sub)
|
|
|
|
# put those data in global 3d array
|
|
data_glob[offset_r:offset_r+nr_sub, offset_t:offset_t+nt_sub, offset_p:offset_p+np_sub] = data_sub
|
|
grid_glob_p[offset_r:offset_r+nr_sub, offset_t:offset_t+nt_sub, offset_p:offset_p+np_sub] = grid_sub_p
|
|
grid_glob_t[offset_r:offset_r+nr_sub, offset_t:offset_t+nt_sub, offset_p:offset_p+np_sub] = grid_sub_t
|
|
grid_glob_r[offset_r:offset_r+nr_sub, offset_t:offset_t+nt_sub, offset_p:offset_p+np_sub] = grid_sub_r
|
|
|
|
# update offset
|
|
offset += n_points_total_sub
|
|
|
|
|
|
return data_glob, grid_glob_r, grid_glob_t, grid_glob_p
|
|
|
|
|
|
|
|
|
|
|
|
|
|
if __name__ == '__main__':
|
|
|
|
# examples
|
|
fpath = './OUTPUT_FILES/out_data_sim_0.h5'
|
|
fpath_grid = './OUPUT_FILES/out_data_grid.h5'
|
|
nr = 10
|
|
nt = 10
|
|
np = 10
|
|
ndiv_r = 2
|
|
ndiv_t = 2
|
|
ndiv_p = 2
|
|
Ks_tomoatt = get_data_from_h5(fpath, fpath_grid, "Data/Ks_inv_0000", nr, nt, np, ndiv_r, ndiv_t, ndiv_p, verbose=True) |