import numpy as np import math from .helper import position_from_rank, rank_from_position, eval_procslice, key_grid from .field import Field3d 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 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 delete(self,key): '''Removes a field from memory.''' import gc 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=None): assert axis<3, "'axis' must be one of 0,1,2." if key_out is None: 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) def qcriterion(self,key_out='Q',from_pressure=False,dtype=np.double): '''Computes Q-criterion''' 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,dtype=dtype) keyx = ('x','y','z') keyu = ('u','v','w') for ii in range(3): key_der = 'd'+keyu[ii]+'d'+keyx[ii] self.differentiate(keyu[ii],ii,key_der) self.shift_to_pressure_grid(key_der) self.field[key_out] += self.field[key_der]**2 self.delete(key_der) # iip1 = (ii+1)%3 key_der1 = 'd'+keyu[ii]+'d'+keyx[iip1] key_der2 = 'd'+keyu[iip1]+'d'+keyx[ii] self.differentiate(keyu[ii],iip1,key_der1) self.differentiate(keyu[iip1],ii,key_der2) self.shift_to_pressure_grid(key_der1) self.shift_to_pressure_grid(key_der2) self.field[key_out] += 2.0*self.field[key_der1]*self.field[key_der2] self.delete(key_der1) self.delete(key_der2) self.field[key_out] *= -0.5 return def shift_to_pressure_grid(self,key,key_out=None): from .field import Field3d if key_out is None: if key in ('u','v','w'): raise ValueError("In-place shift of {'u','v','w'} is not permitted.") key_out = key fld_ref = Field3d.pseudo_field(self.dim,self.origin,self.spacing) 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) 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])