cached stencils

examples/advection_gpu.py

@@ -139,7 +139,7 @@ if __name__=='__main__':
                 'use_builtin_copy':True,
                 'stretching': {
                     'formulation': StretchingFormulation.GRAD_UW,
-                    'order':2
+                    'order':64
                 }
             }
         }

hysop/codegen/maths/stencil/stencil.py

-import math, hashlib
+import os, math, hashlib, gzip, pickle
 import numpy as np
 import scipy as sp
 import sympy as sm
@@ -9,9 +9,34 @@ import scipy.sparse
 import gmpy2 as gmp
 from gmpy2 import mpq
 
+from hysop.tools.io_utils import IO
+
 def _mpqize(x):
     return mpq(x.p,x.q)
 mpqize = np.vectorize(_mpqize)
+    
+_hysop_cache_folder = IO.default_cache_path()
+_stencil_cache_file = _hysop_cache_folder + '/' + 'stencils.pklz'
+
+def _load_cached_stencils():
+    if not os.path.exists(_hysop_cache_folder):
+        os.makedirs(_hysop_cache_folder)
+    if not os.path.isfile(_stencil_cache_file):
+        return {}
+    else:
+        with gzip.open(_stencil_cache_file, 'rb') as f:
+            return pickle.load(f)
+    
+def _export_stencils():
+    with gzip.open(_stencil_cache_file, 'wb') as f:
+        pickle.dump(_cached_stencils,f)
+
+def _clear_stencil_cache():
+    _cached_stencils = {}
+    if os.path.isfile(_stencil_cache_file):
+        os.remove(_stencil_cache_file)
+
+_cached_stencils = _load_cached_stencils()
 
 class Stencil(object):
     def __init__(self, stencil_data, origin, order):
@@ -76,23 +101,6 @@ class Stencil(object):
         return not np.any(self.stencil*mask)
 
 
-
-class StencilGenerator(object):
-    def __init__(self,derivative,dim,dtype=mpq):
-        self.derivative = derivative
-        self.dim        = dim
-        self.dtype      = dtype
-
-        
-        #L = L if np.isscalar(L) else np.asarray(L))
-        #if R is None:
-            #R = k+order-L-1
-            #assert((R>=0).all(), 'R<0...')
-        #else:
-            #assert((L+R+1 == k+order).all()), 'L+R+1 != k+order, cannot compute coefficients !'
-            #R = R if np.isscalar(R) else np.asarray(R))
-        #N = L+R+1
-
 class Stencil1D(Stencil):
     def __init__(self):
         self.Lx    = 0
@@ -101,10 +109,16 @@ class Stencil1D(Stencil):
         self.order = 0
         self.N     = 0
         self.dtype = None
+
+    @staticmethod
+    def _hash_stencil_key(Lx,Rx,k,order,dtype):
+        s = 'Stencil1D_{}_{}_{}_{}_{}'.format(Lx,Rx,k,order,dtype)
+        return hashlib.sha256(s).hexdigest()
+
     
     # generate a discrete stencil of the 'k'-th derivative at order 'order' defined over [-Lx,Rx] grid points. 
     # Note: generated order can be > specified order if stencil coefficients are 0 at borders.
-    def generate(self,Lx,Rx,k,order,dtype=mpq):
+    def generate(self,Lx,Rx,k,order,dtype=mpq,override_cache=False):
         assert(Lx+Rx+1 == k+order), 'Lx+Rx != k+order, cannot compute coefficients !'
         self.Lx    = Lx
         self.Rx    = Rx
@@ -113,19 +127,31 @@ class Stencil1D(Stencil):
         self.dtype = dtype
 
         N = Lx+Rx+1
-        A = np.zeros((N,N),dtype=object)
-        b = np.zeros(N,dtype=object)
-        for i in xrange(N):
-            b[i] = dtype(i==k)
-            for j in xrange(0,N):
-                A[i,j] = dtype(j-Lx)**i
-        A = sm.Matrix(N,N,A.flatten())
-        b = sm.Matrix(N,1,b.flatten())
-
-        S = A.LUsolve(b)
-        S *= dtype(math.factorial(k))
-        S = mpqize(np.asarray(S.tolist(),dtype=object).ravel())
         
+        skey = Stencil1D._hash_stencil_key(Lx,Rx,k,order,dtype)
+        if '1D' not in _cached_stencils:
+            _cached_stencils['1D'] = {}
+
+        if (not override_cache) and (skey in _cached_stencils['1D']):
+            S = _cached_stencils['1D'][skey]
+        else:
+            A = np.zeros((N,N),dtype=object)
+            b = np.zeros(N,dtype=object)
+            for i in xrange(N):
+                b[i] = dtype(i==k)
+                for j in xrange(0,N):
+                    A[i,j] = dtype(j-Lx)**i
+            A = sm.Matrix(N,N,A.flatten())
+            b = sm.Matrix(N,1,b.flatten())
+
+            S = A.LUsolve(b)
+            S *= dtype(math.factorial(k))
+            S = mpqize(np.asarray(S.tolist(),dtype=object).ravel())
+
+            _cached_stencils['1D'][skey] = S
+            _export_stencils()
+
+            
         self.N       = N
         self.stencil = S
         return self