From b91f5e95bb4fb7751a054afd527628d58a872325 Mon Sep 17 00:00:00 2001 From: Michael Krayer Date: Fri, 25 Mar 2022 16:39:14 +0100 Subject: [PATCH] further improved serial --- field.py | 19 ++++++++- serial.py | 121 ++++++++++++++++++++++++++++++++++++++++++------------ 2 files changed, 112 insertions(+), 28 deletions(-) diff --git a/field.py b/field.py index d4464f0..4d83889 100644 --- a/field.py +++ b/field.py @@ -202,7 +202,7 @@ class Field3d: r._dim = dim return r - def save(self,file,name='Field3d',truncate=False): + def save(self,file,name='Field3d',truncate=False,single_precision=False): import h5py is_open = isinstance(file,(h5py.File,h5py.Group)) flag = 'w' if truncate else 'a' @@ -210,7 +210,10 @@ class Field3d: g = f.create_group(name) g.create_dataset('origin',data=self.origin) g.create_dataset('spacing',data=self.spacing) - g.create_dataset('data',data=self.data) + if single_precision: + g.create_dataset('data',data=self.data,dtype='f4') + else: + g.create_dataset('data',data=self.data) if not is_open: f.close() return @@ -590,6 +593,18 @@ class Field3d: val += c[ii,jj,kk]*self.data[i0+ii,j0+jj,k0+kk] return val + def scipy_interpolate(self,pts,oob_fill=None): + from scipy.interpolate import RegularGridInterpolator + if oob_fill is None: + bounds_error = True + fill_value = None + else: + bounds_error = False + fill_value = oob_fill + rgi = RegularGridInterpolator(self.grid(),self.data,method='linear', + bounds_error=bounds_error,fill_value=fill_value) + return rgi(pts) + @staticmethod def padding(origin,spacing,dim,padding,numpad): if isinstance(numpad,int): diff --git a/serial.py b/serial.py index 6687f9b..6365f53 100644 --- a/serial.py +++ b/serial.py @@ -2,6 +2,7 @@ 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""" @@ -17,6 +18,7 @@ class SPP: self.spacing = tuple(spacing) self.dim = tuple(dim) self.xperiodic,self.yperiodic,self.zperiodic = periodicity + self.__KEYTMP = '__TMP' return @classmethod @@ -90,19 +92,39 @@ class SPP: 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=None): + 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 = 'd'+key+('dx','dy','dz')[axis] + 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' @@ -114,9 +136,15 @@ class SPP: 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,dtype=np.double): + 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) @@ -124,40 +152,82 @@ class SPP: 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) + 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] - 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) + 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] - 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) + 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 - return + 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: - if key in ('u','v','w'): - raise ValueError("In-place shift of {'u','v','w'} is not permitted.") - key_out = key + 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) - rel_shift = self.field[key].relative_shift(fld_ref) - self.field[key_out] = self.field[key].shift_origin(rel_shift) + 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 @@ -199,5 +269,4 @@ class SPP: sl_out[axis] = idx_ self.field[key].data[tuple(sl_out)] = self.field[key].data[tuple(sl_in)] idx_ += 1 - return - + return \ No newline at end of file