further improved serial

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):

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:
+        if key_out is None: key_out = self.__KEYTMP
-            key_out = 'd'+key+('dx','dy','dz')[axis]
+        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,
+            self.field[key_out] = Field3d.allocate(self.dim,self.origin,self.spacing,fill=0.0)
-                                                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)
+                if not key_der in self.field:
-                self.shift_to_pressure_grid(key_der)
+                    self.differentiate(keyu[ii],ii,key_out=key_der)
-                self.field[key_out] += self.field[key_der]**2
+                self.shift_to_pressure_grid(key_der,key_out='der')
-                self.delete(key_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)
+                if not key_der1 in self.field:
-                self.differentiate(keyu[iip1],ii,key_der2)
+                    self.differentiate(keyu[ii],iip1,key_out=key_der1)
-                self.shift_to_pressure_grid(key_der1)
+                self.shift_to_pressure_grid(key_der1,key_out='der1')
-                self.shift_to_pressure_grid(key_der2)
+                if not keep_derivatives: self.delete(key_der1)
-                self.field[key_out] += 2.0*self.field[key_der1]*self.field[key_der2]
+                if not key_der2 in self.field:
-                self.delete(key_der1)
+                    self.differentiate(keyu[iip1],ii,key_out=key_der2)
-                self.delete(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_out is None: key_out = self.__KEYTMP
-            if key in ('u','v','w'):
+        if key_out==key:
-                raise ValueError("In-place shift of {'u','v','w'} is not permitted.")
+            raise ValueError("In-place shift 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)
+        if self.field[key].has_same_grid(fld_ref):
-        self.field[key_out] = self.field[key].shift_origin(rel_shift)
+            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
@@ -200,4 +270,3 @@ class SPP:
                 self.field[key].data[tuple(sl_out)] = self.field[key].data[tuple(sl_in)]
                 idx_ += 1
         return
-