From d6acd372e1b33daddfc857bf578174869e9c60ec Mon Sep 17 00:00:00 2001
From: Jean-Baptiste Keck <Jean-Baptiste.Keck@imag.fr>
Date: Fri, 5 Apr 2019 17:19:14 +0200
Subject: [PATCH] polynomial
---
hysop/numerics/interpolation/interpolation.py | 13 +-
hysop/numerics/interpolation/polynomial.py | 243 ++++++++++++++++++
2 files changed, 254 insertions(+), 2 deletions(-)
create mode 100644 hysop/numerics/interpolation/polynomial.py
diff --git a/hysop/numerics/interpolation/interpolation.py b/hysop/numerics/interpolation/interpolation.py
index 0a5bd47e1..be6dc9a6c 100644
--- a/hysop/numerics/interpolation/interpolation.py
+++ b/hysop/numerics/interpolation/interpolation.py
@@ -2,8 +2,17 @@
from hysop.tools.enum import EnumFactory
Interpolation = EnumFactory.create('Interpolation',
- ['LINEAR', 'CUBIC', 'CHEBYSHEV',
- 'L4_4', 'M4', 'Mp4'])
+ ['LINEAR', # requires 0 ghosts
+ 'CUBIC_FDC2', # requires 1 ghosts (estimate derivatives with 2nd order centered fd)
+ 'CUBIC_FDC4', # requires 2 ghosts (estimate derivatives with 4th order centered fd)
+ 'QUINTIC_FDC2', # requires 2 ghosts
+ 'QUINTIC_FDC4' # requires 3 ghosts
+ 'SEPTIC_FDC2', # requires 3 ghosts
+ 'SEPTIC_FDC4' # requires 4 ghosts
+ 'CUBIC', # derivatives order is specified with SpaceDiscretization instead
+ 'QUINTIC', # derivatives order is specified with SpaceDiscretization instead
+ 'SEPTIC', # derivatives order is specified with SpaceDiscretization instead
+ 'L4_4', 'M4', 'Mp4']) # scales interpolators
class MultiScaleInterpolation(object):
pass
diff --git a/hysop/numerics/interpolation/polynomial.py b/hysop/numerics/interpolation/polynomial.py
new file mode 100644
index 000000000..71c97caa0
--- /dev/null
+++ b/hysop/numerics/interpolation/polynomial.py
@@ -0,0 +1,243 @@
+import itertools as it
+import numpy as np
+import scipy as sp
+import sympy as sm
+#import sympy.abc, sympy.polys, sympy.solvers, sympy.functions
+
+from hysop.numerics.stencil.stencil_generator import StencilGenerator, MPQ
+from hysop.tools.sympy_utils import tensor_xreplace
+
+class PolynomialInterpolator(object):
+
+ def __init__(self, dim, deg, fd_order, verbose=False):
+ """
+ Create a PolynomialInterpolator.
+
+ Parameters
+ ----------
+ dim: int
+ Number of dimensions to interpolate.
+ deg: int
+ Polynomial degree (1=linear, 3=cubic, 5=quintic, 7=septic, ...)
+ Degree should be odd.
+ fd_order: int
+ Order of centered finite differences used to compute derivatives.
+ Will affect the number of ghosts of the method.
+ Should be even because this interpolator only use centered dinite differences.
+ verbose: bool
+ Enable or disabled verbosity, default to False.
+
+ Attributes
+ ----------
+ dim: int
+ Number of dimensions to interpolate.
+ deg: int
+ Polynomial degree (1=linear, 3=cubic, 5=quintic, 7=septic, ...)
+ Degree should be odd: deg=2k+1
+ fd: int
+ Order of centered finite differences used to compute derivatives.
+ p: int
+ Corresponds to deg+1.
+ The number of polynomial coefficients corresponds to p**dim.
+ k: int
+ Max derivative order required to compute the polynomial interpolator coefficients.
+ Corresponds to (deg-1)/2.
+ linear: 0
+ cubic: 1
+ quintic: 2
+ septic: 3
+ nonic: 4
+ ghosts: int
+ Return the number of ghosts required by the interpolator on each axis.
+ Corresponds to (k>0)*[fd_order//2 - 1 + (k+1)//2]
+ deg k (k+1)/2 | FDC2 FDC4 FDC6
+ linear: 1 0 0 | 0 0 0
+ cubic: 3 1 1 | 1 2 3
+ quintic: 5 2 1 | 1 2 3
+ septic: 7 3 2 | 2 3 4
+ nonic: 9 4 2 | 2 3 4
+ n: int
+ Corresponds to 2*ghosts+1, the number of required nodes to generate the
+ polynomial coefficients (in each direction).
+
+ See Also
+ --------
+ :class:`PolynomialPatchInterpolator`: Precompute weights for fixed subgrid interpolation.
+ """
+
+ assert dim>0, 'dim<=0'
+ assert deg%2==1, 'deg % 2 != 1'
+ assert fd_order%2==0, 'fd_order % 2 != 0'
+
+ dim = dim
+ k = (deg-1)/2
+
+ if (k>0):
+ required_ghosts = (k+1)/2
+ ghosts = (fd_order/2) - 1 + required_ghosts
+ else:
+ ghosts = 0
+
+ p = deg+1
+ n = 2*ghosts+1
+
+ self.dim = dim
+ self.deg = deg
+ self.p = p
+ self.k = k
+ self.n = n
+ self.fd = fd_order
+ self.ghosts = ghosts
+ self.verbose = verbose
+
+ self._collect_stencils()
+ self._build_system()
+
+ def _collect_stencils(self):
+ dim = self.dim
+ ghosts = self.ghosts
+ verbose = self.verbose
+ order = self.fd
+
+ k = self.k
+ n = self.n
+ ghosts = self.ghosts
+
+ if verbose:
+ print '\nCollecting 1D stencils:'
+
+ S = []
+ SG = StencilGenerator()
+ for i in xrange(k+1):
+ if (k==0):
+ Si = (0,)*ghosts + (1,) + (0,)*ghosts
+ else:
+ Si = SG.generate_exact_stencil(dx=1, dim=1, dtype=MPQ,
+ derivative=i, origin=ghosts, order=order, shape=n)
+ print Si
+ S.append(Si.stencil)
+ if verbose:
+ print ' {}-th derivative: {}'.format(i,Si)
+ self.S = S
+
+ def _build_stencil(self, dvec):
+ dvec = np.asarray(dvec)
+ assert dvec.size == self.dim, 'dvec.size != dim'
+ assert (dvec>=0).all(), 'dvec < 0 => {}'.format(dvec)
+ assert (dvec<=self.k).all(), 'dvec > dmax => {}'.format(dvec)
+
+ Sd = self.S[dvec[0]].copy()
+ for di in dvec[1:]:
+ Sdi = self.S[di]
+ Sd = np.tensordot(Sd, Sdi, axes=0)
+ for i in xrange(0,self.dim-1):
+ Sd = np.flip(Sd, axis=i)
+ return Sd
+
+ def _build_system(self):
+ dim = self.dim
+ deg = self.deg
+ maxd = self.k
+ ghosts = self.ghosts
+ verbose = self.verbose
+
+ forigin = (ghosts,)*dim
+ fshape = (2*(ghosts+1),)*dim
+ pshape = (deg+1,)*dim
+
+ xvals, xvars = tensor_symbol('x',shape=(dim,))
+ fvals, fvars = tensor_symbol('F',fshape,forigin)
+ pvals, pvars = tensor_symbol('C',pshape)
+
+ self.xvals, self.xvars = xvals, xvars
+ self.fvals, self.fvars = fvals, fvars
+ self.pvals, self.pvars = pvals, pvars
+
+ P0 = 0
+ for idx in it.product(range(0,deg+1), repeat=dim):
+ val = pvals[idx]
+ for i in xrange(dim):
+ val*= pow(xvals[i],idx[i])
+ P0 += val
+ self.P0 = P0
+
+ if verbose:
+ print '\nGenerating variables:'
+ print ' *space vars: '
+ print xvals
+ print ' *grid values:'
+ print fvals
+ print ' *polynomial coefficients:'
+ print pvals
+ print ' *polynomial patch:'
+ print P0
+ print '\nBuilding system...'
+
+ eqs = []
+ for dvec in it.product( range(0,maxd+1), repeat=dim ):
+ if verbose:
+ print ' => derivative {}'.format(dvec)
+
+ dP0 = P0
+ for i,deg in enumerate(dvec):
+ dP0 = sm.diff(dP0, xvals[i], deg)
+
+ stencil = self._build_stencil(dvec)
+
+ for idx in it.product( range(ghosts,ghosts+2), repeat=dim):
+ if verbose:
+ print ' -> point {}'.format(idx)
+
+ pos = np.asarray(idx)-ghosts
+ pos = dict(zip(xvals, pos))
+ eq = dP0.xreplace(pos)
+
+ for offset in it.product( range(-ghosts, ghosts+1), repeat=dim):
+ fidx = tuple(np.asarray(idx)+np.asarray(offset))
+ sidx = tuple(np.asarray(offset)+ghosts)
+ eq -= fvals[fidx]*stencil[sidx]
+
+ eqs.append(eq)
+ if verbose:
+ print ' {}'.format(eq)
+
+ if verbose:
+ print '\nSolving system...'
+
+ sol = sm.solve(eqs, pvars, simplify=False)
+ for pvar in pvars:
+ psol = sol[pvar]
+ if verbose:
+ print '{} = {}'.format(pvar,psol)
+ self.sol = sol
+
+ if verbose:
+ print '\nBuilding matrix...'
+ M = np.empty(shape=(pvals.size,fvals.size), dtype=object)
+ for i,pvar in enumerate(pvars):
+ psol = sol[pvar]
+ for j,fvar in enumerate(fvars):
+ M[i,j] = psol.coeff(fvar)
+ self.M = M
+ self.Mflt = M.astype(float)
+ if verbose:
+ original = np.get_printoptions()
+ np.set_printoptions(precision=2, suppress=True, linewidth=1000, threshold=10000)
+ np.set_printoptions(**original)
+ #print self.Mflt
+
+ def interpolate(self, fvals):
+ fvals = np.asarray(fvals)
+ assert fvals.shape == self.fvals.shape
+ #fvars = dict(zip(self.fvars,fvals.ravel()))
+ #cvars = self.sol.copy()
+ #for k,v in cvars.iteritems():
+ #cvars[k] = v.xreplace(fvars)
+ pvals = self.Mflt.dot(fvals.ravel())
+ P0 = self.P0.xreplace(dict(zip(self.pvars, pvals)))
+ return sm.utilities.lambdify(self.xvars[::-1], P0)
+
+
+
+if __name__ == '__main__':
+ P1 = PolynomialInterpolator(dim=2, deg=1, fd_order=2, verbose=True)
--
GitLab