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+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])