diff --git a/hysop/codegen/maths/__init__.py b/hysop/codegen/maths/__init__.py
deleted file mode 100644
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000
diff --git a/hysop/codegen/maths/interpolation/__init__.py b/hysop/codegen/maths/interpolation/__init__.py
deleted file mode 100644
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000
diff --git a/hysop/codegen/maths/interpolation/kernel_generator.py b/hysop/codegen/maths/interpolation/kernel_generator.py
deleted file mode 100644
index 2509599ef59edb7c262e649e4c042d66cfe9872a..0000000000000000000000000000000000000000
--- a/hysop/codegen/maths/interpolation/kernel_generator.py
+++ /dev/null
@@ -1,386 +0,0 @@
-
-import stencil
-import os, hashlib, pickle
-
-import numpy as np
-import matplotlib.pyplot as plt
-
-import scipy as sp
-import scipy.interpolate
-
-import sympy as sm
-import sympy.abc, sympy.polys, sympy.solvers, sympy.functions
-
-import gmpy2 as gmp
-from gmpy2 import mpq,mpfr,f2q
-mpqize = np.vectorize(lambda x: f2q(x))
-
-
-user_home = os.path.expanduser('~')
-hysop_cache_folder = user_home + '/.cache/hysop'
-kernels_cache_file = hysop_cache_folder + '/' + 'kernels.bin'
-
-def _load_cached_kernels():
- if not os.path.exists(hysop_cache_folder):
- os.makedirs(hysop_cache_folder)
- if not os.path.isfile(kernels_cache_file):
- return {}
- else:
- with open(kernels_cache_file, 'rb') as f:
- return pickle.load(f)
-
-def _export_kernels():
- with open(kernels_cache_file, 'wb') as f:
- pickle.dump(_cached_kernels,f)
-
-def _clear_kernel_cache():
- _cached_kernels = {}
- if os.path.isfile(kernels_cache_file):
- os.remove(kernels_cache_file)
-
-def _hash_kernel_key(n,r,deg,Ms,H,remesh):
- s = '{}_{}_{}_{}_{}_{}'.format(n,r,deg,Ms,H,int(remesh))
- return hashlib.sha256(s).hexdigest()
-
-_cached_kernels = _load_cached_kernels()
-
-
-class Kernel(object):
-
- def __init__(self,register=False,verbose=False,split_polys=False,**kargs):
-
- dic = kargs
- for var in ['n','r','deg','Ms','Mh','H','remesh','P']:
- setattr(self,var,dic[var])
-
- if register:
- self._register(dic)
-
- self._build(verbose,split_polys)
-
- def _build(self,verbose,split_polys):
-
- #polynom symbolic variables
- x = sm.abc.x
- t = sm.abc.t
-
- #usefull vars
- Ms = self.Ms
- deg = self.deg
- P = self.P
-
- if verbose:
- print '\t=> Substitution in Polynomials'
- for Pix in P:
- print '\t',Pix.all_coeffs()[::-1]
- print
- for Pix in P:
- print '\t',sm.horner(Pix)
- print
-
- #split polynomials
- X = np.arange(-Ms,+Ms+1)
- if split_polys:
- Pt_l = []
- Pt_r = []
- Cl = np.empty(shape=(deg+1,2*Ms),dtype=float)
- Cr = np.empty(shape=(deg+1,2*Ms),dtype=float)
- for i,Pix in enumerate(P):
- Pit_l = Pix.xreplace({x:t+X[i]})
- Pit_r = Pix.xreplace({x:X[i]+1-t})
- Cl[:,i] = np.asarray(Pit_l.all_coeffs(), dtype=float)
- Cr[:,i] = np.asarray(Pit_r.all_coeffs(), dtype=float)
- Pt_l.append(Pit_l)
- Pt_r.append(Pit_r)
-
- if verbose:
- print '\t=> Splitting polynomials'
- for p in Pt_l:
- print '\t',p.all_coeffs()
- print
- for p in Pt_r:
- print '\t',p.all_coeffs()
- print
- print
- for p in Pt_l:
- print '\t',sm.horner(p)
- print
- for p in Pt_r:
- print '\t',sm.horner(p)
- else:
- C = np.empty(shape=(deg+1,2*Ms),dtype=float)
- for i,Pix in enumerate(P):
- C[:,i] = np.asarray(Pix.all_coeffs(), dtype=float)
-
- if verbose:
- print
- print '\t=> Generating lambdas'
-
- if split_polys:
- gamma_l = sp.interpolate.PPoly(Cl,X,extrapolate=False)
- gamma_r = sp.interpolate.PPoly(Cr,X,extrapolate=False)
- def gamma(x):
- i = np.floor(x)
- z = x-i
- return (z>=0.5)*gamma_r(i+1-z) + (z<0.5)*gamma_l(x)
- else:
- gamma = sp.interpolate.PPoly(C,X,extrapolate=False)
-
- self.I = X
- self.gamma = gamma
- if split_polys:
- self.gamma_l = gamma_l
- self.gamma_r = gamma_r
- self.poly_splitted = True
- else:
- self.poly_splitted = False
-
- if verbose:
- print '\tAll done!'
-
- def _register(self,dic):
- h = _hash_kernel_key(self.n,self.r,self.deg,self.Ms,self.H,self.remesh)
- if h in _cached_kernels.keys():
- raise KeyError('Current kernel was already cached or kernels key clash...')
- _cached_kernels[h] = dic
- _export_kernels()
-
- def __call__(self,x):
- Ms = self.Ms
- return np.logical_and(x>-Ms,x<Ms)*self.gamma(x)
-
-
-class KernelGenerator(object):
-
- def __init__(self,verbose=False):
- self.verbose = verbose
- self.configured = False
-
- def configure(self,n,H=None):
- assert n>0 and n%2==0, 'n not even or n<=0!'
- Ms = n//2+1
-
- if H is None:
- H = np.zeros(2*Ms+1,dtype=int)
- H[Ms] = 1
- self.remesh = True
- else:
- H = mpqize(np.asarray(H))
- assert H.size == 2*Ms+1
- self.remesh = False
- self.H = H
-
- if self.remesh:
- self.Mh = None
- else:
- Mh = []
- for q in xrange(n+1):
- mq = mpq(0)
- for i,h in enumerate(H,start=-Ms):
- mq += (i**q) * h
- Mh.append(mq)
- self.Mh = Mh
-
- self.n = n
- self.Ms = Ms
- self.configured = True
- return self
-
-
- def solve(self, r, override_cache=False):
- assert self.configured, 'KernelGenerator has not been configured!'
- assert r>=0, 'r<0!'
- deg = 2*r+1
-
- n = self.n
- Ms = self.Ms
- H = self.H
- Mh = self.Mh
-
- remesh = self.remesh
- verbose = self.verbose
-
-
- if verbose:
- print 'KernerlGenerator::Solve'
- print '\tn = {}'.format(n)
- print '\tr = {}'.format(r)
- print '\tMs = {}'.format(Ms)
- print '\tdeg = {}'.format(deg)
- print
- print '\tH = ['+','.join([h.__str__() for h in H])+']'
- if not self.remesh:
- print '\tMh = ['+','.join([m.__str__() for m in Mh])+']'
- print
-
- H = mpqize(np.asarray(H))
- if not self.remesh:
- Mh = mpqize(np.asarray(Mh))
-
- # check if kernel was not already cached
- if override_cache:
- if verbose:
- print '\tCache overwrite requested!'
- else:
- if verbose:
- print '\tChecking if kernel was already computed... ',
- h = _hash_kernel_key(n,r,deg,Ms,H,remesh)
- if h in _cached_kernels:
- if verbose:
- print 'True'
- print '\tLoading kernel from cache!'
- print
- kernel = Kernel(verbose=verbose,register=False,**_cached_kernels[h])
- return kernel
- elif verbose:
- print 'False'
- print '\tBuilding linear system!'
- print
-
- #polynom symbolic variable
- x = sm.abc.x
-
- #build Ms*(deg+1) symbolic coefficients (polynomial unknowns)
- coeffs = []
- for k in xrange(Ms):
- coeffs.append([sm.symbols('C{k}_{d}'.format(k=k,d=d)) for d in xrange(deg+1)])
- #build discrete moments rhs values
- M = []
- for i in xrange(n+1):
- Mi = []
- for j in xrange(deg+1):
- if remesh and i==j:
- Mij = f2q(1)
- elif not remesh and i-j>=0 and Mh[i-j]!=0:
- Mij = sm.symbols('M{i}_{j}'.format(i=i,j=j))
- else:
- Mij = 0
- Mi.append(Mij)
- M.append(Mi)
-
- #build the Ms polynomials
- Pp = []
- Pm = []
- for i,C in enumerate(coeffs):
- pexpr = f2q(0)
- mexpr = f2q(0)
- for d,c in enumerate(C):
- pexpr += c*(x**d)
- mexpr += c*((-x)**d)
- Pp.append(sm.polys.polytools.poly(pexpr,x))
- Pm.append(sm.polys.polytools.poly(mexpr,x))
- P = Pm[::-1] + Pp
- Pbeg = -Ms
- Pend = +Ms-1
- Prange = lambda: xrange(-Ms,+Ms)
- Penum = lambda: enumerate(P,start=Pbeg)
-
- #precompute the r first polynomial derivatives
- dPs = []
- for p in Pp:
- dP = [p]
- for i in xrange(r):
- p = p.diff()
- dP.append(p)
- dPs.append(dP)
-
- # Build overdetermined system of equations
- eqs = []
-
- ### Continuity equations (Gamma is Cr continuous)
- # Parity in x=0 -- Gamma is EVEN -> (r+1)//2 eqs
- for d in xrange(1,r+1,2):
- eq = coeffs[0][d] # =0
- eqs.append(eq)
- # Right-most point, zero derivatives -> (r+1) eqs
- for d in xrange(r+1):
- eq = dPs[-1][d].xreplace({x:Ms}) # =0
- eqs.append(eq.as_expr())
- # Cr-continuity on inner points -> (Ms-1)*(r+1) eqs
- for d in xrange(r+1):
- for i in xrange(Ms-1):
- eq = dPs[i][d].xreplace({x:i+1}) - dPs[i+1][d].xreplace({x:i+1}) # = 0
- eqs.append(eq.as_expr())
-
- ### Interpolation condition on the left -> Ms equations
- for i in xrange(Ms):
- eq = Pp[i].xreplace({x:i}) - H[Ms+i] # = 0
- eqs.append(eq.as_expr())
-
- ### Discrete moments
- s = sm.symbols('s')
- for m in xrange(0,n+1):
- expr = f2q(0)
- for l in xrange(-Ms+1,Ms+1):
- if m>0 and l==0: continue
- i = Ms-l
- e = P[i].xreplace({x:s-f2q(l)})
- if m>0: e *= f2q(l**m)
- expr += e
- Pm = sm.polys.polytools.poly(expr,s)
- for i,Cmi in enumerate(Pm.all_coeffs()[::-1]):
- eqs.append(Cmi-M[m][i])
-
- if verbose:
- print
- print '\t=> System built.'
-
- unknowns = [c for cl in coeffs for c in cl]+[m for ml in M for m in ml if isinstance(m,sm.Symbol)]
- if verbose:
- print '\tUnknowns: ',unknowns
-
- sol = sm.solve(eqs,unknowns)
- if len(sol)!=len(unknowns):
- if verbose:
- print 'sol=',sol
- raise RuntimeError('Could not solve linear system!')
- elif verbose:
- print
- print '\t=> System solved.'
- print '\t',sol
- print
-
- for i,Pix in enumerate(P):
- P[i] = Pix.xreplace(sol)
-
- kernel = Kernel(n=n,r=r,deg=deg,Ms=Ms,
- H=H, Mh=Mh, remesh=remesh,
- P=P, register=True, verbose=verbose)
- return kernel
-
-
-if __name__=='__main__':
- verbose=True
-
- #H = stencil.SymmetricStencil1D().generate(L=3,k=2,order=1).add_zeros(1,1)
- for i in xrange(1,5):
- p = 2*i
- kg = KernelGenerator(verbose).configure(p)
- kernels = []
- for r in [1,2,4]:#,8,16,32]:
- try:
- kernels.append(kg.solve(r))
- except:
- pass
-
- if len(kernels)==0:
- continue
- continue
- k0 = kernels[0]
-
- fig = plt.figure()
- plt.xlabel(r'$x$')
- plt.ylabel(r'$\Lambda_{'+'{},{}'.format(p,'r')+'}$')
- X = np.linspace(-k0.Ms,+k0.Ms,1000)
- s = plt.subplot(1,1,1)
- for i,k in enumerate(kernels):
- s.plot(X,k(X),label=r'$\Lambda_{'+'{},{}'.format(p,k.r)+'}$')
- s.plot(k0.I,k0.H,'or')
- axe_scaling = 0.10
- ylim = s.get_ylim()
- Ly = ylim[1] - ylim[0]
- s.set_ylim(ylim[0]-axe_scaling*Ly,ylim[1]+axe_scaling*Ly)
- s.legend()
-
- plt.show()
- plt.savefig('out/gamma_{}_r'.format(p))
diff --git a/hysop/codegen/maths/runge_kutta/__init__.py b/hysop/codegen/maths/runge_kutta/__init__.py
deleted file mode 100644
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000
diff --git a/hysop/codegen/maths/runge_kutta/runge_kutta.py b/hysop/codegen/maths/runge_kutta/runge_kutta.py
deleted file mode 100644
index 8f01365bfbc030e14f78207b39a20ebf44aaf37b..0000000000000000000000000000000000000000
--- a/hysop/codegen/maths/runge_kutta/runge_kutta.py
+++ /dev/null
@@ -1,111 +0,0 @@
-
-import numpy as np
-from hysop.codegen.runge_kutta_coeffs import Itype,Qtype,rk_params, available_methods
-
-
-# default methods to dump rational and integers expression to string
-def Q_dump(q):
- return '({}/{})'.format(q.numerator,q.denominator)
-def I_dump(i):
- return '{}'.format(z)
-
-
-class RungeKutta(object):
- pass
-
-# Tni = Tn + alpha_i*dt
-# Xni = Xn + dt*sum(j=0,i-1,gamma_ij*Knj)
-# Kni = F(Tni,Xni)
-# X_n+1 = Xn + dt*sum(i=0,M-1,beta_i*Ki)
-class ExplicitRungeKutta(RungeKutta):
-
- def __init__(self,method,I_dump=I_dump,Q_dump=Q_dump):
-
- if method not in available_methods():
- raise ValueError('{} was not implemented yet! Valid values are {}.'.format(name,implemented_methods.keys()))
-
- params = rk_params[method]
-
- self.method = method
-
- self.order = params['order']
- self.stages = params['stages']
-
- self.alpha = params['alpha']
- self.beta = params['beta']
- self.gamma = params['gamma']
-
- self.I_dump = I_dump
- self.Q_dump = Q_dump
-
- def dump(self,val):
- if isinstance(val,Itype):
- return self.I_dump(val)
- elif isinstance(val,Qtype):
- return self.Q_dump(val)
- else:
- return '{}'.format(val)
-
- # Tni = Tn + alpha_i*dt
- def Tni(self,i,Tn,dt):
- alpha = self.alpha[i]
- if alpha == 0:
- return '{}'.format(Tn)
- elif alpha == 1:
- return '{} + {}'.format(Tn,dt)
- else:
- return '{} + {}*{}'.format(Tn,self.dump(alpha),dt)
-
- # Xni = Xn + dt*sum(j=0,i-1,gamma_ij*Knj)
- def Xni_sum(self,i):
- _sum = ''
- gamma = self.gamma[i-1,:]
- for j in xrange(i):
- g = gamma[j]
- if g == 0:
- continue
- elif g == 1:
- _sum += 'K[{}]'.format(j)
- else:
- _sum += '{}*K[{}]'.format(self.dump(g),j)
- _sum += ' + '
- _sum = _sum[:-3]
- return _sum
- def Xni(self, i,Xn,dt):
- if i>0:
- _sum = Xni_sum(i)
- return '{} + {}*({})'.format(Xn,dt,_sum)
- else:
- return '{}'.format(Xn)
-
- # X_n+1 = Xn + dt*sum(i=0,M-1,beta_i*Ki)
- def step_sum(self):
- _sum = ''
- for i in xrange(self.stages):
- beta = self.beta[i]
- if beta == 0:
- continue
- elif beta == 1:
- _sum += 'K[{}]'.format(i)
- else:
- _sum += '{}*K[{}]'.format(self.dump(beta),i)
- _sum += ' + '
- _sum = _sum[:-3]
- return _sum
- def step(self,Xn,dt):
- return '{} + {}*({})'.format(Xn,dt,self.step_sum())
-
-
-
-if __name__ == '__main__':
- R = ExplicitRungeKutta('RK4')
- for i in xrange(4):
- print R.Tni(i,Tn='To',dt='dt')
- print
- for i in xrange(4):
- print R.Xni(i,Xn='Xo',dt='dt')
- print
- print R.step(Xn='Xo',dt='dt')
-
-
-
diff --git a/hysop/codegen/maths/runge_kutta/runge_kutta_coeffs.py b/hysop/codegen/maths/runge_kutta/runge_kutta_coeffs.py
deleted file mode 100644
index 3f952fc0fda212bb764fb4836ac5de7fb981ffe0..0000000000000000000000000000000000000000
--- a/hysop/codegen/maths/runge_kutta/runge_kutta_coeffs.py
+++ /dev/null
@@ -1,70 +0,0 @@
-
-import numpy as np
-import gmpy2 as gmp
-I = gmp.mpz
-Q = gmp.mpq
-Itype = I(0).__class__
-Qtype = Q(1,2).__class__
-
-Z = I(0)
-
-# For all i in {0,...,M-1}
-# Tni = Tn + alpha_i*dt
-# Xni = Xn + dt*sum(j=0,i-1,gamma_ij*Kj)
-# Kni = F(Tni,Xni)
-# X_n+1 = Xn + dt*sum(i=0,M-1,beta_i*Ki)
-rk_params = {}
-
-RK1 = {}
-RK1['order'] = 1
-RK1['stages'] = 1
-RK1['alpha'] = [Z]
-RK1['beta'] = [I(1)]
-RK1['gamma'] = [[]]
-rk_params['RK1'] = RK1
-rk_params['euler'] = RK1
-
-RK2 = {}
-RK2['order'] = 2
-RK2['stages'] = 2
-RK2['alpha'] = [Z,Q(1,2)]
-RK2['beta'] = [Z,I(1)]
-RK2['gamma'] = [[Q(1,2)]]
-rk_params['RK2'] = RK2
-
-
-RK4 = {}
-RK4['order'] = 4
-RK4['stages'] = 4
-RK4['alpha'] = [Z,Q(1,2),Q(1,2),I(1)]
-RK4['beta'] = [Q(1,6), Q(1,3), Q(1,3), Q(1,6)]
-RK4['gamma'] = [[Q(1,2), Z, Z ],
- [Z, Q(1,2), Z ],
- [Z, Z, I(1)]]
-rk_params['RK4'] = RK4
-
-
-# check and wrap in numpy arrays
-for m in rk_params:
- C = rk_params[m]
- S = C['stages']
-
- alpha = np.asarray(C['alpha'], dtype=object)
- beta = np.asarray(C['beta'], dtype=object)
- gamma = np.asarray(C['gamma'], dtype=object)
-
- assert alpha.shape == (S,)
- assert beta.shape == (S,)
- if S>1:
- assert gamma.shape == (S-1,S-1,)
- else:
- assert gamma.size == 0
-
- C['alpha'] = alpha
- C['beta'] = beta
- C['gamma'] = gamma
-
- rk_params[m] = C
-
-def available_methods():
- return rk_params.keys()
diff --git a/hysop/codegen/maths/stencil/__init__.py b/hysop/codegen/maths/stencil/__init__.py
deleted file mode 100644
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000
diff --git a/hysop/codegen/maths/stencil/coeffs.py b/hysop/codegen/maths/stencil/coeffs.py
deleted file mode 100644
index d670935400d437de476ce9e2eeefebdccc580437..0000000000000000000000000000000000000000
--- a/hysop/codegen/maths/stencil/coeffs.py
+++ /dev/null
@@ -1,9 +0,0 @@
-
-import numpy as np
-import gmpy2 as gmp
-I = gmp.mpz
-Q = gmp.mpq
-Itype = I(0).__class__
-Qtype = Q(1,2).__class__
-
-Z = I(0)
diff --git a/hysop/codegen/maths/stencil/solver.py.todo b/hysop/codegen/maths/stencil/solver.py.todo
deleted file mode 100644
index 7c61125290d82fed9252f3c86d33653eb773685d..0000000000000000000000000000000000000000
--- a/hysop/codegen/maths/stencil/solver.py.todo
+++ /dev/null
@@ -1,275 +0,0 @@
-
-
-import itertools as it
-import operator, math, sys
-
-import numpy as np
-import sympy as sm
-from sympy.core.singleton import S as singleton
-from sympy.solvers.inequalities import reduce_abs_inequality, solve_poly_inequality, reduce_rational_inequalities, reduce_inequalities
-
-x,y = sm.symbols('x,y')
-f = sm.symbols('f')
-
-dec = ['\u208{}'.format(i).decode('unicode-escape') for i in xrange(10)]
-sign = [u'\u208a',u'\u208b']
-parenthesis = [u'\u208d', u'\u208e']
-
-def subscript(i,with_sign=False):
- decimals = '0123456789'
- snumber = str(i)
- if with_sign:
- s0 = snumber[0]
- if s0 in decimals:
- snumber = '+'+snumber
- out = u''
- for s in snumber:
- if s in decimals:
- out += dec[int(s)]
- elif s=='+':
- out += sign[0]
- elif s=='-':
- out += sign[1]
- else:
- out += s
- return out
-
-def subscripts(ids,sep,with_sign=False,with_parenthesis=False,prefix=''):
- if with_parenthesis:
- return u'{}{}{}{}'.format(prefix,parenthesis[0],sep.join([subscript(i,with_sign) for i in ids]),parenthesis[1])
- else:
- return u'{}{}'.format(prefix,sep.join([subscript(i,with_sign) for i in ids]))
-
-def tensor_symbol(prefix,shape,origin=None,sep=None,with_parenthesis=False,force_sign=False):
- origin = origin if origin is not None else np.asarray([0]*len(shape))
- sep = sep if sep is not None else u','
-
- with_sign = force_sign or ((origin>0).any() and len(shape)>1)
- tensor = np.empty(shape=shape,dtype=object)
- for idx in np.ndindex(*shape):
- ids = idx-origin
- sname = subscripts(ids,sep,with_sign=with_sign,with_parenthesis=with_parenthesis,prefix=prefix)
- tensor[idx] = sm.Symbol(sname.encode('utf-8'), real=True)
- tensor_vars = tensor.ravel().tolist()
- return tensor, tensor_vars
-
-def tensor_xreplace(tensor,vars):
- T = tensor.copy()
- for idx in np.ndindex(*tensor.shape):
- symbol = tensor[idx]
- if (symbol in vars.keys()):
- T[idx] = vars[symbol]
- elif (symbol.name in vars.keys()):
- T[idx] = vars[symbol.name]
- return T
-
-def factor_split(expr,variables,constant_var=None,include_var=False,init=0,_factor=True,_handle_const=True):
- expr = expr.expand()
- factors = {}
- for var in variables:
- factors[var] = init
- for arg in expr.args:
- I = arg.atoms(sm.Symbol).intersection(variables)
- if len(I)==0:
- continue
- elif len(I)==1:
- var = I.pop()
- if not include_var:
- arg = arg.xreplace({var:1})
- factors[var] += arg
- else:
- assert False, 'Expression containing two or more variables!\n{} contains {}.'.format(arg,I)
- #for var in factors:
- #factors[var] = sm.factor(factors[var],var)
- #cterm = expr
- #for var in factors:
- #cterm -= factors[var]
- #if constant_var is None:
- #assert cterm==0, 'No const var provided but there is a constant term \'{}\'!'.format(cterm)
- #else:
- #factors[constant_var] = m
- #if not include_var:
- #for var in factors:
- #factors[var] = sm.subs(factors[var],{var:1})
- return factors
-
-def build_eqs_from_dicts(d0,d1):
- treated = []
- eqs = []
- for k in d0.keys():
- expr0 = d0[k]
- expr1 = d1[k] if k in d1.keys() else 0
- de = expr1-expr0
- if de != 0:
- eqs.append(de)
- treated.append(k)
- for k in d1.keys():
- if k not in treated:
- expr0 = 0
- expr1 = d1[k]
- de = expr1-expr0
- if de != 0:
- eqs.append(de)
- return eqs
-
-def prod(it,init=1):
- return reduce(operator.mul, it, init)
-
-def gen_vec(basename,dim,same=False):
- if same:
- vec = np.array([sm.Symbol(basename, real=True)])
- else:
- vec,_ = tensor_symbol(prefix=basename,shape=(dim,))
- return vec
-
-
-def taylor(df,dx,maxn):
- expr = 0
- for n in xrange(maxn+1):
- expr += taylorn(df,dx,n)
- return expr
-
-def taylorn(df,dx,n):
- dim = dx.size
-
- def preficate(it):
- return sum(it)==n
-
- nn = range(n+1)
- itd = it.product(nn,repeat=dim)
- itd = it.ifilter(preficate,itd)
- expr = 0
- for der in itd:
- expr += taylorn_term(df,dx,der)
- return expr
-
-def taylorn_term(df,dx,der):
- fact = np.asarray([math.factorial(d) for d in der])
- return df[der]*prod((dx**der)/fact)
-
-
-
-dim = 3
-order = 1
-L=np.asarray([-2]*dim)
-R=np.asarray([+2]*dim)
-
-origin=np.abs(L)
-window = R-L+1
-dx = gen_vec('dx',dim,same=True)
-
-df,df_vars = tensor_symbol(prefix='f',shape=[9]*dim)
-S,svars = tensor_symbol(prefix='S',shape=window,origin=origin)
-
-expr = 0
-last = 0
-last2 = 0
-last3 = 0
-for ids in np.ndindex(*window):
- offset = ids-origin
- expr += S[ids]*taylor(df,offset*dx,5)
- last += S[ids]*taylorn(df,offset*dx,6)
- last2 += S[ids]*taylorn(df,offset*dx,7)
- last3 += S[ids]*taylorn(df,offset*dx,8)
-fact = factor_split(expr,df_vars)
-user_eqs = {}
-#user_eqs[df[2,0]] = 1
-#user_eqs[df[0,2]] = 1
-user_eqs[df[2,0,0]] = 1
-user_eqs[df[0,2,0]] = 1
-user_eqs[df[0,0,2]] = 1
-
-eqs = build_eqs_from_dicts(fact,user_eqs)
-
-sol = sm.solve(eqs,svars)
-for s in sol:
- print s,' = ',sol[s]
-
-Ss = tensor_xreplace(S,sol)
-
-#sol2 = {}
-#sol2[S[2,2]] = sm.Integer(0)/2/dx[0]**2
-
-
-print 'MORE MORE SOL'
-new_eqs = []
-nvars = set()
-s = 0
-for k,v in factor_split(last,df_vars).iteritems():
- if v!=0:
- eq = sm.simplify(v.xreplace(sol))
- symbols = eq.atoms(sm.Symbol)
- if eq != 0:
- print k, ' -> ', eq
- s+=eq
- inter = symbols.intersection(svars)
- if inter:
- new_eqs.append(eq)
- nvars = nvars.union(inter)
-sol2 = sm.solve(new_eqs,nvars,ignore=dx.tolist())
-for k,v in factor_split(last2,df_vars).iteritems():
- if v!=0:
- eq = sm.simplify(v.xreplace(sol).xreplace(sol2))
- symbols = eq.atoms(sm.Symbol)
- if eq != 0:
- print k, ' -> ', eq
- s+=eq
- inter = symbols.intersection(svars)
- if inter:
- new_eqs.append(eq)
- nvars = nvars.union(inter)
-sol3 = sm.solve(new_eqs,nvars,ignore=dx.tolist())
-print 'SOL:'
-for k,v in sol3.iteritems():
- print k,'=',v
-print 'MORE MORE MORE SOL'
-new_eqs = []
-nvars = set()
-for k,v in factor_split(last3,df_vars).iteritems():
- if v!=0:
- eq = sm.simplify(v.xreplace(sol).xreplace(sol2).xreplace(sol3))
- symbols = eq.atoms(sm.Symbol)
- if eq != 0:
- print k, ' -> ', eq
- inter = symbols.intersection(svars)
- if inter:
- new_eqs.append(eq)
- nvars = nvars.union(inter)
-sol4 = sm.solve(new_eqs,nvars,ignore=dx.tolist())
-print 'SOL:'
-for k,v in sol4.iteritems():
- print k,'=',v
-for s in [sol2,sol3,sol4]:
- for k,e in sol.iteritems():
- sol[k] = e.xreplace(s)
- sol.update(s)
-print sol
-Ss = tensor_xreplace(S,sol)*dx[0]**2
-print Ss
-
-E = np.zeros_like(S)
-k = gen_vec('k',dim)
-cfl = sm.Symbol('c')
-for ids in np.ndindex(*window):
- offset = ids-origin
- E[ids] = sm.exp(sm.I*np.dot(k,offset))
-stability = (1 + cfl*np.sum(Ss*E))
-stability = stability.rewrite(sm.cos)
-
-cosit = it.product([0,1],repeat=dim)
-cond = True
-econd = []
-for cit in cosit:
- subs = {}
- for i in xrange(dim):
- subs[k[i]] = cit[i]*sm.pi
- stab = stability.xreplace(subs).xreplace({dx[0]:1})
- if cfl in stab.atoms(sm.Symbol):
- try:
- cond &= reduce_abs_inequality(sm.Abs(stab)-1,'<',cfl)
- except:
- econd.append(sm.Abs(stab)<1)
-
-print cond.as_set()
-print econd
-
diff --git a/hysop/codegen/maths/stencil/stencil.py b/hysop/codegen/maths/stencil/stencil.py
deleted file mode 100644
index 8205d2fe81a72ab4516aa2c68e9fa18b13deddd4..0000000000000000000000000000000000000000
--- a/hysop/codegen/maths/stencil/stencil.py
+++ /dev/null
@@ -1,168 +0,0 @@
-import math, hashlib
-import numpy as np
-import scipy as sp
-import sympy as sm
-import itertools as it
-
-import scipy.sparse
-
-import gmpy2 as gmp
-from gmpy2 import mpq
-
-def _mpqize(x):
- return mpq(x.p,x.q)
-mpqize = np.vectorize(_mpqize)
-
-class Stencil(object):
- def __init__(self, stencil_data, origin, order):
- stencil = np.asarray(stencil_data)
- origin = np.asarray(origin)
- order = np.asarray(order)
-
- dim = stencil.ndim
- dtype = stencil.dtype
- shape = stencil.shape
-
- assert (origin>=0).all(), 'origin<0!'
- L = origin
- R = shape-origin-1
-
- self.stencil = stencil
- self.origin = origin
- self.order = order
-
- self.dim = dim
- self.dtype = dtype
- self.shape = shape
-
- self.L = L
- self.R = R
-
- def coo_stencil(self):
- assert self.dim==2
- return sp.sparse.coo_matrix(self.stencil,shape=self.shape,dtype=self.dtype)
-
- def iteritems(self):
- iterator = np.ndindex(self.shape)
- iterator = it.imap(lambda x: (x-self.origin,self.stencil[x]), iterator)
- iterator = it.ifilter(lambda x: x[1]!=0, iterator)
- return iterator
-
- def hash(self):
- stencil = self.stencil
- m = hashlib.sha1()
- m.update(self.origin)
- m.update(self.order)
- for d in stencil:
- m.update(str(d))
- return m.hexdigest()[:8]
-
- def is_centered(self,axe=None):
- if axe==None:
- return (L==R).all()
- else:
- return (L[axe]==R[axe]).all()
-
- def is_cross(self):
- mask = np.ones(self.shape,dtype=bool)
- access = self.origin.tolist()
- for i in xrange(self.dim):
- acc = [x for x in access]
- acc[i] = slice(0,self.shape[i])
- mask[acc] = False
- 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
- self.Rx = 0
- self.k = 0
- self.order = 0
- self.N = 0
- self.dtype = None
-
- # 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):
- assert(Lx+Rx+1 == k+order), 'Lx+Rx != k+order, cannot compute coefficients !'
- self.Lx = Lx
- self.Rx = Rx
- self.k = k
- self.order = order
- 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())
-
- self.N = N
- self.stencil = S
- return self
-
- def add_zeros(self,L,R):
- S = self.stencil
- Zl = np.zeros(L,dtype=object)
- Zr = np.zeros(R,dtype=object)
- return np.concatenate((Zl,S,Zr))
-
- def __str__(self):
- if(self.N == 0):
- return ''
- else:
- return '['+','.join([s.__str__() for s in self.stencil])+']'
-
-class SymmetricStencil1D(Stencil1D):
- def __init__(self):
- super(self.__class__,self).__init__()
-
- def generate(self,L,k,order,dtype=mpq):
- assert L == k+order, 'L != k+order, cannot compute stencil!'
- Lx = (L-1)//2
- Ly = Lx
- return super(self.__class__,self).generate(Lx,Ly,k,order,dtype)
-
-if __name__ == '__main__':
- print 'first derivative'
- print Stencil1D().generate(Lx=0,Rx=1,k=1,order=1,dtype=mpq)
- print Stencil1D().generate(Lx=1,Rx=0,k=1,order=1,dtype=mpq)
- print Stencil1D().generate(Lx=1,Rx=1,k=1,order=2,dtype=mpq)
- print Stencil1D().generate(Lx=2,Rx=2,k=1,order=4,dtype=mpq)
- print
- print 'second derivative'
- print Stencil1D().generate(Lx=1,Rx=1,k=2,order=1,dtype=mpq)
- print Stencil1D().generate(Lx=1,Rx=2,k=2,order=2,dtype=mpq)
- print Stencil1D().generate(Lx=2,Rx=1,k=2,order=2,dtype=mpq)
- print Stencil1D().generate(Lx=2,Rx=2,k=2,order=3,dtype=mpq)
- print
- print 'using symmetric stencil class'
- for i in xrange(1,5):
- print SymmetricStencil1D().generate(L=2*i+1,k=2,order=2*(i-1)+1)
-
diff --git a/hysop/codegen/maths/sympy_utils.py b/hysop/codegen/maths/sympy_utils.py
deleted file mode 100644
index 8e7cd3604f4b9cfbf0b9c245a3c6308765a83334..0000000000000000000000000000000000000000
--- a/hysop/codegen/maths/sympy_utils.py
+++ /dev/null
@@ -1,45 +0,0 @@
-
-import sympy as sm
-
-def non_eval_xreplace(expr, rule):
- """
- Duplicate of sympy's xreplace but with non-evaluate statement included
- """
- if expr in rule:
- return rule[expr]
- elif rule:
- args = []
- altered = False
- for a in expr.args:
- try:
- new_a = non_eval_xreplace(a, rule)
- except AttributeError:
- new_a = a
- if new_a != a:
- altered = True
- args.append(new_a)
- args = tuple(args)
- if altered:
- return expr.func(*args, evaluate=False)
- return expr
-
-# Convert powers to mult. in polynomial expressions
-# Example: x^3 -> x*x*x
-def remove_pows(expr):
- pows = list(expr.atoms(sm.Pow))
- repl = [sm.Mul(*[b]*e,evaluate=False) for b,e in [i.as_base_exp() for i in pows]]
- e = non_eval_xreplace(expr,dict(zip(pows,repl)))
- print e
- return e
-
-def evalf_str(x,n,literal='',significant=True):
- x = x.evalf(n).__str__()
- if significant:
- i=len(x)
- while i>1 and x[i-1] == '0':
- i-=1
- if i>1:
- x = x[:i+1]
- return x+literal
-
-
diff --git a/hysop/numerics/finite_differences.py b/hysop/numerics/finite_differences.py
index 42bc646cabede1d0d018d1b515cd35fdef9af2c4..6ae41654b90d271590590e5c57d63f03bf0610c6 100644
--- a/hysop/numerics/finite_differences.py
+++ b/hysop/numerics/finite_differences.py
@@ -48,6 +48,7 @@ Note
"""
from abc import ABCMeta, abstractmethod
from hysop.constants import debug
+from hysop.numerics.stencil import Stencil, StencilGenerator
import hysop.tools.numpywrappers as npw
import numpy as np
diff --git a/hysop/numerics/stencil.py b/hysop/numerics/stencil.py
index cf179c9a07166910d632dc02835d6f8d85808578..ceb9655095992e951871796992ccf3629483da8c 100644
--- a/hysop/numerics/stencil.py
+++ b/hysop/numerics/stencil.py
@@ -1,3 +1,9 @@
+"""Finite differences stencil inteface used in various numerical methods.
+
+* :class:`~hysop.numerics.stencil.Stencil`
+* :class:`~hysop.numerics.stencil.StencilGenerator`
+
+"""
import math, hashlib
import itertools as it
@@ -6,23 +12,32 @@ import numpy as np
import sympy as sm
import scipy as sp
-import scipy.sparse
+import scipy.linalg
import gmpy2 as gmp
from gmpy2 import mpq
+
+MPQ = mpq().__class__
def _mpqize(x):
- return mpq(x.p,x.q)
+ x = mpq(x.p,x.q)
+ return x
mpqize = np.vectorize(_mpqize)
+
class Stencil(object):
"""
- Stencil used for finite abitrary dimension differences.
+ Stencil used for finite abitrary dimension differences.
+ """
- Parameters
+ def __init__(self, coeffs, origin, order):
+ """
+ Stencil used for finite abitrary dimension differences.
+
+ Attributes
----------
- coeffs: array_like
+ coeffs: np.ndarray
Coefficients of the stencil.
- origin: array_like
+ origin: np.ndarray
Origin of the stencil from left-most point.
order: int
Order of the stencil.
@@ -36,11 +51,6 @@ class Stencil(object):
Distance from origin to left-most point included in stencil support.
R: np.ndarray
Distance from origin to right-most point included in stencil support.
- """
-
- def __init__(self, coeffs, origin, order):
- """
- Stencil used for finite abitrary dimension differences.
Parameters
----------
@@ -54,13 +64,14 @@ class Stencil(object):
Raises
------
ValueError
- Raised if one of component of `origin` is negative.
+ Raised if one component of `origin` is negative.
See Also
--------
- StencilGenerator: Generate Stencil objects
- """
- stencil = np.asarray(stencil_data)
+ :class:`StencilGenerator`: Generate Stencil objects.
+ """
+
+ stencil = np.asarray(coeffs)
origin = np.asarray(origin)
order = np.asarray(order)
@@ -87,16 +98,30 @@ class Stencil(object):
def coo_stencil(self):
"""
Return a 2d stencil as a sparse coordinate matrix.
- @return : scipy.sparse.coo_matrix
+
+ Returns
+ -------
+ coo_stencil: scipy.sparse.coo_matrix
+ Sparse Coordinate Matrix of the stencil.
+
+ Raises
+ ------
+ RuntimeError
+ Raised if stencil is not 2d.
"""
- assert self.dim==2
+ if self.dim!=2:
+ raise RuntimeError('Stencil is not 2d !')
return sp.sparse.coo_matrix(self.stencil,shape=self.shape,dtype=self.dtype)
def iteritems(self):
"""
- Return an (offset,coefficient) iterator iterating on all non zero coefficients.
+ Return an (offset,coefficient) iterator iterating on all **non zero** coefficients.
Offset is taken from origin.
- @return : itertools.iterator
+
+ Returns
+ -------
+ it : itertools.iterator
+ Zipped offset and coefficient iterator.
"""
iterator = np.ndindex(self.shape)
iterator = it.imap(lambda x: (x-self.origin,self.stencil[x]), iterator)
@@ -108,8 +133,16 @@ class Stencil(object):
Hash the stencil with its origin, order and coefficients.
The hash algorithm used is sha1.
By default only the 8 first chars are kept from the generated hash.
- @parameter keep_only : number of chars kept from the hashed hexedecimal string.
- @return : hexadecimal hash string
+
+ Parameters
+ ----------
+ keep_only : int
+ Number of chars kept from the hashed hexedecimal string.
+
+ Returns
+ -------
+ hash : string
+ Hexadecimal hash string of size :attr:`keep_only`.
"""
stencil = self.stencil
m = hashlib.sha1()
@@ -122,7 +155,10 @@ class Stencil(object):
def is_centered(self,axe=None):
"""
Check if the stencil is centered.
- @return : bool
+
+ Returns
+ -------
+ is_centered : bool
"""
if axe==None:
return (L==R).all()
@@ -132,7 +168,10 @@ class Stencil(object):
def is_cross(self):
"""
Check if the stencil is a cross (zeros everywhere except on a nd cross).
- @return : bool
+
+ Returns
+ -------
+ is_cross : bool
"""
mask = np.ones(self.shape,dtype=bool)
access = self.origin.tolist()
@@ -142,26 +181,130 @@ class Stencil(object):
mask[acc] = False
return not np.any(self.stencil*mask)
+ def __str__(self):
+ return '''
+== {}D Stencil ==
+origin: {}
+order: {}
+shape: {}
+dtype: {}
+data: {}
+=================
+'''.format(self.dim,self.origin,self.order,self.shape,self.dtype,self.stencil)
+
+
class StencilGenerator(object):
"""
- Generate a stencil from its order and derivative order.
- Results are cached and compressed locally.
+ Generate various stencils.
+ Results are cached and compressed locally.
"""
- def __init__(self,dim,dtype=mpq):
+ def __init__(self, dim, dtype=MPQ, cache_override=False):
"""
- Initialize a StencilGenerator to generate stencils of dimension :attr:`dim` and type :attr:`dtype`
- Results are cached and compressed locally.
+ StencilGenerator to generate stencils of dimension :attr:`dim` and type :attr:`dtype`.
+ Results are cached and compressed locally.
+
+ Attributes
+ ----------
+ dim : int
+ Dimension of the generated stencils.
+ dtype : {np.float32, np.float64, mpq}
+ Datatype of the generated stencils.
+
+ Parameters
+ ----------
+ dim : int
+ Dimension of the generated stencils.
+ dtype : {np.float32, np.float64, hysop.numerics.stencil.MPQ}
+ Datatype of the generated stencils.
+
+ Other Parameters
+ ----------------
+ cache_override : bool
+ If set, overwrite already cached stencils.
+
+ Raises
+ ------
+ ValueError
+ Raised if dim<1 or if dtype in [np.float32, np.float64] and dim!=1.
+ TypeError
+ Raised if dtype is not a valid type.
+
+ Notes
+ -----
+ Depending on data type :attr:`dtype` a different solver is used.
+ For fast generation of 1D stencils, use `np.float32` or `np.float64`.
+ For n-dimentional exact fractional coefficients stencils use `mpq` (default).
+ Exact fractional results are cached.
+
+ See Also
+ --------
+ :class:`~Stencil`
"""
- self.derivative = derivative
+
+ if dim<1:
+ raise ValueError('dim<1!')
+ if dtype not in [np.float32, np.float64, MPQ]:
+ raise TypeError('Invalid type!')
+ if dtype in [np.float32, np.float64] and dim!=1:
+ raise ValueError('Bad configuration!')
+
self.dim = dim
self.dtype = dtype
- L = np.asarray([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
+ def generate(self, origin, derivative, min_order):
+ dim = self.dim
+ min_order = np.asarray([min_order]*dim, dtype=int) if np.isscalar(min_order) else np.asarray(min_order)
+ origin = np.asarray([origin]*dim, dtype=int) if np.isscalar(origin) else np.asarray(origin)
+ derivative = np.asarray([derivative]*dim, dtype=int) if np.isscalar(derivative) else np.asarray(derivative)
+ if (min_order.size != dim) or (origin.size != dim) or (derivative.size != dim):
+ raise ValueError('Bad input dimensions!')
+ if (min_order<0).any() or (origin<0).any() or (derivative<0).any():
+ raise ValueError('Bad input values!')
+
+ if self.dtype in [np.float32, np.float64]:
+ stencil = self._generate_approximative_stencil(origin,derivative,min_order)
+ else:
+ stencil = self._generate_exact_stencil(origin,derivative,min_order)
+ return stencil
+
+ def _generate_approximative_stencil(self, origin, derivative, min_order):
+ dim = self.dim
+ dtype = self.dtype
+ assert dim==1
+
+ origin = origin[0]
+ derivative = derivative[0]
+ min_order = min_order[0]
+
+ k = derivative
+ Lx = origin
+ Rx = k+min_order-Lx-1
+ while Rx<0:
+ min_order+=1
+ Rx = k+min_order-Lx-1
+
+ N = Lx+Rx+1
+ A = np.zeros((N,N),dtype=dtype)
+ b = np.zeros(N,dtype=dtype)
+ for i in xrange(N):
+ b[i] = dtype(i==k)
+ for j in xrange(N):
+ A[i,j] = dtype(j-Lx)**i
+
+ try:
+ S = scipy.linalg.solve(A,b,overwrite_a=True,overwrite_b=True)
+ S *= dtype(math.factorial(k))
+ return Stencil(S,origin,min_order)
+ except:
+ print 'Cannot generate stencil (singular system), trying with the exact solver!'
+ return self._generate_exact_stencil(origin,derivative,min_order)
+
+ def _generate_exact_stencil(self, origin, derivative, min_order):
+ return None
+
+if __name__ == '__main__':
+ S0 = StencilGenerator(dim=1,dtype=np.float64)
+ for k in xrange(1,100):
+ print S0.generate(k,1,2*k)
+