272 lines
12 KiB
Python
272 lines
12 KiB
Python
import numpy as np
|
|
import math
|
|
from .helper import position_from_rank, rank_from_position, eval_procslice, key_grid
|
|
from .field import Field3d
|
|
import warnings
|
|
|
|
class SPP:
|
|
"""Serial Python Postprocessor for suspend"""
|
|
def __init__(self,field,periodicity,bounds,origin,spacing,dim):
|
|
assert all([key in field for key in ('u','v','w')]), "Required fields: 'u','v','w'."
|
|
for axis in range(3):
|
|
if periodicity[axis]:
|
|
assert bounds[2*axis]==0.0, "Lower bound in periodic direction must equal zero."
|
|
self.field = field
|
|
self.bounds = tuple(bounds)
|
|
self.periodicity = tuple(periodicity)
|
|
self.origin = tuple(origin)
|
|
self.spacing = tuple(spacing)
|
|
self.dim = tuple(dim)
|
|
self.xperiodic,self.yperiodic,self.zperiodic = periodicity
|
|
self.__KEYTMP = '__TMP'
|
|
return
|
|
|
|
@classmethod
|
|
def from_snapshot(cls,snap,procslice=None,dtype=np.double,load_pressure=False,verbose=False):
|
|
# Load fields
|
|
keys = ('u','v','w','p') if load_pressure else ('u','v','w')
|
|
field = {}
|
|
for key_ in keys:
|
|
if verbose:
|
|
print('[SPP.from_snapshot] loading {}'.format(key_))
|
|
field[key_] = Field3d.from_snapshot(snap,key_,procslice,dtype)
|
|
# As a common grid, the pressure grid is used. Extract its offset
|
|
proc = snap.procgrid()
|
|
origin = snap.origin()
|
|
spacing = snap.spacing()
|
|
orip,dimp = eval_procslice(proc,origin,spacing,'p',procslice)[1:]
|
|
orip = tuple(orip[ii]-spacing[ii] for ii in range(3))
|
|
dimp = tuple(dimp[ii]+2 for ii in range(3))
|
|
# Set periodicity and global bounds
|
|
bounds = snap.bounds()
|
|
periodicity = list(snap.periodicity())
|
|
if procslice is not None: # if not all procs are loaded, set periodicity to false
|
|
if procslice[0] is not None or procslice[1] is not None:
|
|
periodicity[0] = False
|
|
if procslice[2] is not None or procslice[3] is not None:
|
|
periodicity[1] = False
|
|
if procslice[4] is not None or procslice[5] is not None:
|
|
periodicity[2] = False
|
|
return cls(field,periodicity,bounds,orip,spacing,dimp)
|
|
|
|
@classmethod
|
|
def from_file(cls,file,verbose=False):
|
|
'''Load saved state from file.'''
|
|
import h5py
|
|
if verbose:
|
|
print('[SPP.from_file] file={:s}'.format(file))
|
|
f = h5py.File(file,'r')
|
|
bounds = tuple(f['bounds'])
|
|
periodicity = tuple(f['periodicity'])
|
|
origin = tuple(f['origin'])
|
|
spacing = tuple(f['spacing'])
|
|
dim = tuple(f['dim'])
|
|
# Load fields
|
|
field_names = tuple(x.decode('ascii') for x in f['field_names'])
|
|
field = {}
|
|
for key_ in field_names:
|
|
if verbose: print('[SPP.from_file] loading {:s}'.format(key_))
|
|
field[key_] = Field3d.from_file(f,key_)
|
|
f.close()
|
|
return cls(field,periodicity,bounds,origin,spacing,dim)
|
|
|
|
def save(self,file,verbose=False):
|
|
'''Save current state to file.'''
|
|
import h5py
|
|
if verbose:
|
|
print('[SPP.save] {:s}'.format(file))
|
|
f = h5py.File(file,'w')
|
|
f.create_dataset('bounds',data=self.bounds)
|
|
f.create_dataset('periodicity',data=self.periodicity)
|
|
f.create_dataset('origin',data=self.origin)
|
|
f.create_dataset('spacing',data=self.spacing)
|
|
f.create_dataset('dim',data=self.dim)
|
|
# Save fields and field names
|
|
ascii_type = h5py.string_dtype('ascii',32)
|
|
d = f.create_dataset('field_names',len(self.field),dtype=ascii_type)
|
|
ii_ = 0
|
|
for key_ in self.field:
|
|
self.field[key_].save(f,name=key_)
|
|
d[ii_] = key_
|
|
ii_ += 1
|
|
f.close()
|
|
return
|
|
|
|
def save_fields(self,file,keys,verbose=False,on_pressure_grid=False,single_precision=False):
|
|
'''Save current state to file.'''
|
|
import h5py
|
|
if not isinstance(keys,(list,tuple,np.ndarray)):
|
|
keys = tuple(keys)
|
|
if verbose: print('[SPP.save_field] {:s}'.format(file))
|
|
f = h5py.File(file,'w')
|
|
for key_ in keys:
|
|
if verbose: print('[SPP.save_field] saving {:s}'.format(key_))
|
|
assert key_ in self.field, \
|
|
"Field {} cannot be saved, because it does not exist.".format(key_)
|
|
if on_pressure_grid:
|
|
tmp = self.shift_to_pressure_grid(key_,key_out=None)
|
|
else:
|
|
tmp = self.field[key]
|
|
tmp.save(f,name=key_,single_precision=single_precision)
|
|
f.close()
|
|
return
|
|
|
|
def delete(self,key):
|
|
'''Removes a field from memory.'''
|
|
import gc
|
|
assert not key in ('u','v','w'), "Deleting u,v,w is not permitted."
|
|
assert key in self.field, \
|
|
"Field {} cannot be deleted, because it does not exist.".format(key)
|
|
del self.field[key]
|
|
gc.collect()
|
|
return
|
|
|
|
def differentiate(self,key,axis,key_out=''):
|
|
assert axis<3, "'axis' must be one of 0,1,2."
|
|
if key_out is None: key_out = self.__KEYTMP
|
|
if key_out=='': key_out = 'd'+key+('dx','dy','dz')[axis]
|
|
shifting_state = self.shifting_state(key,axis=axis)
|
|
if shifting_state==-1:
|
|
shift_origin = 'after'
|
|
elif shifting_state==0:
|
|
shift_origin = 'before'
|
|
elif shifting_state==1:
|
|
shift_origin = 'before'
|
|
else:
|
|
raise ValueError("Invalid shifting state.")
|
|
self.field[key_out] = self.field[key].derivative(axis,shift_origin=shift_origin)
|
|
self.exchange_ghost_cells(key_out)
|
|
if key_out==self.__KEYTMP:
|
|
tmp = self.field[key_out]
|
|
self.delete(key_out)
|
|
return tmp
|
|
return self.field[key_out]
|
|
|
|
def qcriterion(self,key_out='Q',from_pressure=False,keep_derivatives=False):
|
|
'''Computes Q-criterion'''
|
|
if key_out is None: key_out = self.__KEYTMP
|
|
if from_pressure:
|
|
assert 'p' in self.field, "Pressure field is not loaded."
|
|
self.field[key_out] = 0.5*self.field['p'].laplacian(only_keep_interior=False)
|
|
else:
|
|
assert 'u' in self.field, "Velocity field is not loaded: u"
|
|
assert 'v' in self.field, "Velocity field is not loaded: v"
|
|
assert 'w' in self.field, "Velocity field is not loaded: w"
|
|
self.field[key_out] = Field3d.allocate(self.dim,self.origin,self.spacing,fill=0.0)
|
|
keyx = ('x','y','z')
|
|
keyu = ('u','v','w')
|
|
for ii in range(3):
|
|
key_der = 'd'+keyu[ii]+'d'+keyx[ii]
|
|
if not key_der in self.field:
|
|
self.differentiate(keyu[ii],ii,key_out=key_der)
|
|
self.shift_to_pressure_grid(key_der,key_out='der')
|
|
self.field[key_out] += self.field['der']**2
|
|
self.delete('der')
|
|
if not keep_derivatives: self.delete(key_der)
|
|
#
|
|
iip1 = (ii+1)%3
|
|
key_der1 = 'd'+keyu[ii]+'d'+keyx[iip1]
|
|
key_der2 = 'd'+keyu[iip1]+'d'+keyx[ii]
|
|
if not key_der1 in self.field:
|
|
self.differentiate(keyu[ii],iip1,key_out=key_der1)
|
|
self.shift_to_pressure_grid(key_der1,key_out='der1')
|
|
if not keep_derivatives: self.delete(key_der1)
|
|
if not key_der2 in self.field:
|
|
self.differentiate(keyu[iip1],ii,key_out=key_der2)
|
|
self.shift_to_pressure_grid(key_der2,key_out='der2')
|
|
if not keep_derivatives: self.delete(key_der2)
|
|
self.field[key_out] += 2.0*self.field['der1']*self.field['der2']
|
|
self.delete('der1')
|
|
self.delete('der2')
|
|
self.field[key_out] *= -0.5
|
|
if key_out==self.__KEYTMP:
|
|
tmp = self.field[key_out]
|
|
self.delete(key_out)
|
|
return tmp
|
|
return self.field[key_out]
|
|
|
|
def vorticity(self,axis,key_out='',keep_derivatives=False):
|
|
keyx = ('x','y','z')
|
|
keyu = ('u','v','w')
|
|
iip1 = (axis+1)%3
|
|
iip2 = (axis+2)%3
|
|
key_der1 = 'd'+keyu[iip1]+'d'+keyx[iip2]
|
|
key_der2 = 'd'+keyu[iip2]+'d'+keyx[iip1]
|
|
if key_out is None: key_out = self.__KEYTMP
|
|
if key_out=='': key_out = 'om'+keyx[axis]
|
|
if not key_der1 in self.field:
|
|
self.differentiate(keyu[iip1],iip2,key_out=key_der1)
|
|
self.shift_to_pressure_grid(key_der1,key_out='der1')
|
|
if not keep_derivatives: self.delete(key_der1)
|
|
if not key_der2 in self.field:
|
|
self.differentiate(keyu[iip2],iip1,key_out=key_der2)
|
|
self.shift_to_pressure_grid(key_der1,key_out='der2')
|
|
if not keep_derivatives: self.delete(key_der2)
|
|
self.field[key_out] = self.field['der2']-self.field['der1']
|
|
self.delete('der1')
|
|
self.delete('der2')
|
|
if key_out==self.__KEYTMP:
|
|
tmp = self.field[key_out]
|
|
self.delete(key_out)
|
|
return tmp
|
|
return self.field[key_out]
|
|
|
|
def shift_to_pressure_grid(self,key,key_out=None):
|
|
from .field import Field3d
|
|
if key_out is None: key_out = self.__KEYTMP
|
|
if key_out==key:
|
|
raise ValueError("In-place shift is not permitted.")
|
|
fld_ref = Field3d.pseudo_field(self.dim,self.origin,self.spacing)
|
|
if self.field[key].has_same_grid(fld_ref):
|
|
self.field[key_out] = self.field[key].copy()
|
|
else:
|
|
rel_shift = self.field[key].relative_shift(fld_ref)
|
|
self.field[key_out] = self.field[key].shift_origin(rel_shift)
|
|
self.exchange_ghost_cells(key_out)
|
|
if key_out==self.__KEYTMP:
|
|
tmp = self.field[key_out]
|
|
self.delete(key_out)
|
|
return tmp
|
|
return self.field[key_out]
|
|
|
|
def shifting_state(self,key,axis=None):
|
|
'''Returns the relative shift with respect to domain
|
|
boundaries as either 1 (= +0.5h), 0 (= 0h) or -1 (= -0.5h).'''
|
|
if axis is None:
|
|
return tuple(self.shifting_state(key,axis=ii) for ii in range(3))
|
|
assert axis<3, "'axis' must be one of 0,1,2."
|
|
return int(round(
|
|
(self.field[key].origin[axis]-self.origin[axis])/(0.5*self.spacing[axis])
|
|
))
|
|
|
|
def exchange_ghost_cells(self,key):
|
|
'''Copy ghost cells in periodic directions.'''
|
|
for axis in range(3):
|
|
if not self.periodicity[axis]:
|
|
continue
|
|
x = self.field[key].grid(axis=axis)
|
|
sl_in = 3*[slice(None)]
|
|
sl_out = 3*[slice(None)]
|
|
# Loop through grid points below lower bound
|
|
xgh = np.mod(x[x<self.bounds[2*axis]],self.bounds[2*axis+1])
|
|
idx_ = 0
|
|
for x_ in xgh:
|
|
ii = self.field[key].nearest_gridpoint(x_,axis=axis)
|
|
assert math.fabs(x_-x[ii])<self.field[key].eps_collapse, \
|
|
"Invalid ghost cell detected: {:f}, {:f}".format(x_,x[ii])
|
|
sl_in[axis] = ii
|
|
sl_out[axis] = idx_
|
|
self.field[key].data[tuple(sl_out)] = self.field[key].data[tuple(sl_in)]
|
|
idx_ += 1
|
|
# Loop through grid points above upper bound
|
|
xgh = np.mod(x[x>=self.bounds[2*axis+1]],self.bounds[2*axis+1])
|
|
idx_ = self.field[key].dim(axis=axis)-len(xgh)
|
|
for x_ in xgh:
|
|
ii = self.field[key].nearest_gridpoint(x_,axis=axis)
|
|
assert math.fabs(x_-x[ii])<self.field[key].eps_collapse, \
|
|
"Invalid ghost cell detected:{:f}, {:f}".format(x_,x[ii])
|
|
sl_in[axis] = ii
|
|
sl_out[axis] = idx_
|
|
self.field[key].data[tuple(sl_out)] = self.field[key].data[tuple(sl_in)]
|
|
idx_ += 1
|
|
return |