diff --git a/examples/multiresolution/scalar_advection.py b/examples/multiresolution/scalar_advection.py
index fa9d3a58601feb4e77c338d357b1040515ed650f..9acfdebc216f35f383f8996061f39952e4fd944b 100644
--- a/examples/multiresolution/scalar_advection.py
+++ b/examples/multiresolution/scalar_advection.py
@@ -26,9 +26,9 @@ def compute(args):
# Define domain
dim = args.ndim
- npts = args.npts
- snpts = args.snpts
- fnpts = tuple(2*_ for _ in snpts)
+ npts = args.npts # coarse resolution
+ snpts = args.snpts # fine resolution
+ fnpts = tuple(3*_ for _ in snpts) # finest resolution
lboundaries = (BoxBoundaryCondition.PERIODIC,)*dim
rboundaries = (BoxBoundaryCondition.PERIODIC,)*dim
box = Box(origin=args.box_origin, length=args.box_length, dim=dim,
diff --git a/hysop/backend/device/opencl/operator/spatial_filtering.py b/hysop/backend/device/opencl/operator/spatial_filtering.py
index c01c11a2b428c596c6d69576ae1c61fbeb2885c9..d6f168e6e38e4f1646a18cbd23ed39ba6d064ce9 100644
--- a/hysop/backend/device/opencl/operator/spatial_filtering.py
+++ b/hysop/backend/device/opencl/operator/spatial_filtering.py
@@ -32,17 +32,17 @@ class OpenClPolynomialInterpolationFilter(PolynomialInterpolationFilterBase, Ope
dim = dFin.dim
assert dFin.is_scalar
assert dFout.is_scalar
+ assert self.subgrid_interpolator.gr == gr
ekg = self.elementwise_kernel_generator
- W = self.subgrid_interpolator.W
- s = self.subgrid_interpolator.s
- n = (self.subgrid_interpolator.n,)*dim
- ghosts = np.asarray((self.subgrid_interpolator.ghosts,)*dim)
+ Wr = self.subgrid_interpolator.Wr
+ n = self.subgrid_interpolator.n
+ ghosts = np.asarray(self.subgrid_interpolator.ghosts)
I = np.asarray(local_indices_symbols[:dim][::-1])
fin, fout = ekg.dfields_to_ndbuffers(dFin, dFout)
Fin = ekg.symbolic_tmp_scalars('F', shape=n, dtype=dFin.dtype)
- Fout_values = W.dot(Fin.ravel()).reshape(s)
+ Fout_values = Wr.dot(Fin.ravel()).reshape(gr)
exprs = ()
for idx in np.ndindex(*n):
@@ -53,7 +53,7 @@ class OpenClPolynomialInterpolationFilter(PolynomialInterpolationFilterBase, Ope
exprs += (e,)
interpolate_grid_kernel, _ = ekg.elementwise_kernel('interpolate_grid',
- *exprs, compute_resolution=dFin.compute_resolution, debug=False)
+ *exprs, compute_resolution=self.iter_shape, debug=False)
exchange_ghosts = self.dFout.exchange_ghosts(build_launcher=True)
diff --git a/hysop/backend/host/python/operator/spatial_filtering.py b/hysop/backend/host/python/operator/spatial_filtering.py
index 24a2fd065f391a26e91b7b47880a2227c3b8cba1..a8268882778ba704c9393adaf7e1c5a0964f49ca 100644
--- a/hysop/backend/host/python/operator/spatial_filtering.py
+++ b/hysop/backend/host/python/operator/spatial_filtering.py
@@ -16,8 +16,7 @@ class PythonPolynomialInterpolationFilter(PolynomialInterpolationFilterBase, Hos
if self.discretized:
return
super(PythonPolynomialInterpolationFilter, self).discretize(**kwds)
- self.W = self.subgrid_interpolator.W.astype(self.dtype)
- self.iter_shape = self.dFin.compute_resolution
+ self.Wr = self.subgrid_interpolator.Wr.astype(self.dtype)
@op_apply
def apply(self, **kwds):
@@ -25,16 +24,17 @@ class PythonPolynomialInterpolationFilter(PolynomialInterpolationFilterBase, Hos
super(PythonPolynomialInterpolationFilter, self).apply(**kwds)
fin = self.fin
fout = self.fout
- shape = self.iter_shape
+ ishape = self.dFin.compute_resolution
+ oshape = self.dFout.compute_resolution
+ periodicity = self.dFin.periodicity
gr = self.grid_ratio
- s, n = self.subgrid_interpolator.s, self.subgrid_interpolator.n
- oshape = fout.shape
+ gr, n = self.subgrid_interpolator.gr, self.subgrid_interpolator.n
+ Wr = self.Wr
- for idx in np.ndindex(*shape):
- oslc = tuple(slice(j*gr[i], (j+1)*gr[i]+1, 1) for i,j in enumerate(idx))
- islc = tuple(slice(1+j, 1+j+n, 1) for i,j in enumerate(idx))
- rslc = tuple(slice(0, min(s[i], oshape[i]-j*gr[i]), 1) for i,j in enumerate(idx))
- fout[oslc] = self.W.dot(fin[islc].ravel()).reshape(s)[rslc]
+ for idx in np.ndindex(*self.iter_shape):
+ oslc = tuple(slice(j*gr[i], (j+1)*gr[i], 1) for i,j in enumerate(idx))
+ islc = tuple(slice(periodicity[i]+j, periodicity[i]+j+n[i], 1) for i,j in enumerate(idx))
+ fout[oslc] = Wr.dot(fin[islc].ravel()).reshape(gr)
self.dFout.exchange_ghosts()
diff --git a/hysop/numerics/interpolation/polynomial.py b/hysop/numerics/interpolation/polynomial.py
index 8aa98499aad8fd9e4da468830edae4a2848f79de..2e49a40aa22714d837920052535e9238568377b7 100644
--- a/hysop/numerics/interpolation/polynomial.py
+++ b/hysop/numerics/interpolation/polynomial.py
@@ -2,6 +2,7 @@ import itertools as it
import numpy as np
import scipy as sp
import sympy as sm
+import warnings
try:
import flint
@@ -16,11 +17,12 @@ except ImportError:
has_flint = False
from hysop.tools.enum import EnumFactory
-from hysop.tools.types import check_instance, InstanceOf
+from hysop.tools.types import check_instance, InstanceOf, to_tuple
from hysop.tools.cache import update_cache, load_data_from_cache
from hysop.tools.io_utils import IO
from hysop.tools.sympy_utils import tensor_xreplace, tensor_symbol
from hysop.tools.decorators import debug
+from hysop.tools.warning import HysopCacheWarning
from hysop.numerics.stencil.stencil_generator import CenteredStencilGenerator, MPQ
from hysop.numerics.interpolation.interpolation import MultiScaleInterpolation
@@ -96,7 +98,7 @@ class PolynomialInterpolationMethod(SpaceDiscretizationMethod):
Operator helper to handle polynomial interpolation method.
"""
__default_method = {
- MultiScaleInterpolation: PolynomialInterpolation.NONIC_FDC2,
+ MultiScaleInterpolation: PolynomialInterpolation.CUBIC_FDC2,
}
__available_methods = {
@@ -168,11 +170,11 @@ class PolynomialInterpolator(object):
----------
dim: int
Number of dimensions to interpolate.
- deg: int
+ deg: int or tuple of ints
Polynomial degree (1=linear, 3=cubic, 5=quintic, 7=septic, ...)
- Degree should be odd.
- fd: int
- Order of centered finite differences used to compute derivatives.
+ Degree should be odd on each axis.
+ fd: int or tuple of ints
+ Order of centered finite differences used to compute derivatives in each direction.
Will affect the number of ghosts of the method.
Should be even because this interpolator only use centered dinite differences.
approximative: bool
@@ -184,20 +186,20 @@ class PolynomialInterpolator(object):
----------
dim: int
Number of dimensions to interpolate.
- deg: int
+ deg: tuple of ints
Polynomial degree (1=linear, 3=cubic, 5=quintic, 7=septic, ...)
Degree should be odd: deg=2k+1
- fd: int
- Order of centered finite differences stencils used to compute derivatives.
- p: int
+ fd: tuple of ints
+ Order of centered finite differences stencils used to compute derivatives for each direction.
+ p: tuple of ints
Corresponds to deg+1.
- The total number of polynomial coefficients corresponds to P=p**dim.
+ The total number of polynomial coefficients corresponds to P=p0*p1*...*p(dim-1).
P: int
- The total number of polynomial coefficients P=p**dim.
- k: int
- Max derivative order required to compute the polynomial interpolator coefficients.
+ The total number of polynomial coefficients P=p0*p1*...*p(dim-1)
+ k: tuple of ints
+ Max derivative order required to compute the polynomial interpolator coefficients in each direction.
Also the regularity of the resulting interpolant. Corresponds to (deg-1)/2.
- ghosts: int
+ ghosts: tuple of ints
Return the number of ghosts required by the interpolator on each axis.
Corresponds to (k>0)*[fd//2 - 1 + (k+1)//2]
deg k (k+1)/2 | FDC2 FDC4 FDC6
@@ -206,29 +208,29 @@ class PolynomialInterpolator(object):
quintic: 5 2 1 | 1 2 3
septic: 7 3 2 | 2 3 4
nonic: 9 4 2 | 2 3 4
- n: int
+ n: tuple of ints
Corresponds to 2*(ghosts+1), the number of required nodes to generate the
polynomial coefficients (in each direction).
- In total we have N=n**dim input nodes.
+ In total we have N=n0*n1*...*n(dim-1) input nodes.
- G G
+ G1 G1
<-> <->
X X X X
X P P X
X P P X
X X X X
<----->
- n
+ n1
N: int
- Total number of input nodes N=n**dim.
+ Total number of input nodes N=n0*n1*...*n(dim-1).
M: np.ndarray
Grid values to polynomial coefficient matrix:
M.dot(F.ravel()) will give C.ravel(), coefficients of P(x0,x1,...)
- N = n**dim
+ N
<--------->
X X X X X X ^ |f0| ^ |c0| ^
X X X X X X | |f1| | |c1| |
- M = X X X X X X | P = p**dim F = |f2| | N C = M*F = |c2| | P
+ M = X X X X X X | P F = |f2| | N C = M*F = |c2| | P
X X X X X X | |f3| | |c3| |
X X X X X X v |f4| | |c4| v
|f5| v
@@ -242,23 +244,44 @@ class PolynomialInterpolator(object):
"""
assert dim>0, 'dim<=0'
- assert deg%2==1, 'deg % 2 != 1'
- assert fd%2==0, 'fd % 2 != 0'
- dim = dim
- k = (deg-1)/2
-
- if (k>0):
- required_ghosts = (k+1)/2
- ghosts = (fd/2) - 1 + required_ghosts
- else:
- ghosts = 0
-
- p = deg+1
- n = 2*(ghosts+1)
+ deg = to_tuple(deg)
+ if len(deg)==1:
+ deg*=dim
+ check_instance(deg, tuple, values=int, size=dim)
- P = p**dim
- N = n**dim
+ fd = to_tuple(fd)
+ if len(fd)==1:
+ fd*=dim
+ check_instance(fd, tuple, values=int, size=dim)
+
+ p = tuple(degi+1 for degi in deg)
+ k = tuple((degi-1)/2 for degi in deg)
+
+ ghosts = ()
+ n = ()
+ for (fdi,ki) in zip(fd,k):
+ if (ki>0):
+ gi = (fdi/2) - 1 + (ki+1)/2
+ else:
+ gi = 0
+ ni = 2*(gi+1)
+ ghosts += (gi,)
+ n += (ni,)
+
+ check_instance(deg, tuple, values=int, size=dim)
+ check_instance(fd, tuple, values=int, size=dim)
+ check_instance(p, tuple, values=int, size=dim)
+ check_instance(k, tuple, values=int, size=dim)
+ check_instance(n, tuple, values=int, size=dim)
+ check_instance(ghosts,tuple, values=int, size=dim)
+ assert all(fdi%2==0 for fdi in fd), 'fd % 2 != 0'
+ assert all(degi%2==1 for degi in deg), 'deg % 2 != 1'
+ assert all(pi%2==0 for pi in p), 'p % 2 != 0'
+ assert all(ni%2==0 for ni in n), 'n % 2 != 0'
+
+ P = np.prod(p, dtype=np.int32)
+ N = np.prod(n, dtype=np.int32)
self.dim = dim
self.deg = deg
@@ -278,7 +301,7 @@ class PolynomialInterpolator(object):
dim = self.dim
ghosts = self.ghosts
verbose = self.verbose
- order = self.fd
+ fd = self.fd
approximative = self.approximative
k = self.k
@@ -289,76 +312,89 @@ class PolynomialInterpolator(object):
print '\nCollecting 1D stencils:'
SG = CenteredStencilGenerator()
- SG.configure(dim=1, dtype=np.float64, order=order)
- S = []
- for i in xrange(k+1):
- msg='Failed to compute stencil derivative={}, order={}, origin={}'
- msg=msg.format(i, order, ghosts)
- try:
- if approximative:
- Si = SG.generate_approximative_stencil(derivative=i)
- else:
- Si = SG.generate_exact_stencil(derivative=i)
- Si.replace_symbols({Si.dx:1})
- except:
- print msg
- raise
- msg+=' got {}.'.format(Si.coeffs)
- assert (not Si.is_symbolic()), msg
- Si = Si.reshape((self.n-1,))
- assert Si.origin == ghosts
- Si = Si.coeffs
- S.append(Si)
+ SG.configure(dim=1, dtype=np.float64)
+ S = {}
+ for direction in xrange(dim):
if verbose:
- print ' {}-th derivative: {}'.format(i,Si)
- self.S = S
+ print ' Direction {}'.format(direction)
+ Sd = S.setdefault(direction, [])
+ nd = n[direction]
+ kd = k[direction]
+ fdd = fd[direction]
+ gd = ghosts[direction]
+ for i in xrange(kd+1):
+ msg='Failed to compute stencil derivative={}, order={}, origin={}'
+ msg=msg.format(i, fdd, gd)
+ try:
+ if approximative:
+ Si = SG.generate_approximative_stencil(order=fdd, derivative=i)
+ else:
+ Si = SG.generate_exact_stencil(order=fdd, derivative=i)
+ Si.replace_symbols({Si.dx:1})
+ except:
+ print msg
+ raise
+ msg+=' got {}.'.format(Si.coeffs)
+ assert (not Si.is_symbolic()), msg
+ Si = Si.reshape((nd-1,))
+ assert Si.origin == gd
+ Si = Si.coeffs
+ Sd.append(Si)
+ if verbose:
+ print ' {}-th derivative: {}'.format(i,Si)
+ return S
def _build_stencil(self, dvec):
dvec = np.asarray(dvec)
+ k = self.k
+ S = self.S
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)
+ assert all(dvec>=0), 'dvec < 0 => {}'.format(dvec)
+ assert all(di<=ki for (di,ki) in zip(dvec,k)), '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)
+ Sd = S[0][dvec[0]].copy()
+ for (d,i) in enumerate(dvec[1:],1):
+ Sd = np.tensordot(Sd, S[d][i], axes=0)
return Sd
def _build_interpolator(self):
- dim = self.dim
- deg = self.deg
- maxd = self.k
+ dim = self.dim
+ deg = self.deg
+ k = self.k
+ n = self.n
+ p = self.p
ghosts = self.ghosts
approximative = self.approximative
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)
+ xvals, xvars = tensor_symbol('x', shape=(dim,))
+ fvals, fvars = tensor_symbol('F', n, ghosts)
+ pvals, pvars = tensor_symbol('C', p)
self.xvals, self.xvars = xvals, xvars
self.fvals, self.fvars = fvals, fvars
self.pvals, self.pvars = pvals, pvars
+ try:
+ data = load_data_from_cache(self.cache_file(), self.key)
+ if (data is not None):
+ (P0, S, M) = data
+ if _check_matrix(M):
+ self.P0 = P0
+ self.S = S
+ self.M = M
+ return
+ except Exception as e:
+ msg='Failed to load data from cache because:\n{}'.format(e)
+ warnings.warn(msg, HysopCacheWarning)
+
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
+ for idx in it.product(*tuple(range(0,pi) for pi in p)):
+ P0 += pvals[idx] * np.prod(np.power(xvals, idx))
self.P0 = P0
- self._collect_stencils()
-
- M = load_data_from_cache(self.cache_file(), self.key)
- if (M is not None) and _check_matrix(M):
- self.M = M
- return
+ S = self._collect_stencils()
+ self.S = S
if verbose:
print '\nGenerating variables:'
@@ -373,7 +409,7 @@ class PolynomialInterpolator(object):
print '\nBuilding system...'
eqs = []
- for dvec in it.product( range(0,maxd+1), repeat=dim ):
+ for dvec in it.product(*tuple(range(0,ki+1) for ki in k)):
if verbose:
print ' => derivative {}'.format(dvec)
@@ -386,17 +422,17 @@ class PolynomialInterpolator(object):
print ' stencil:'
print stencil
- for idx in it.product( range(ghosts,ghosts+2), repeat=dim):
+ for idx in it.product(*tuple(range(gi,gi+2) for gi in ghosts)):
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)
+
+ for offset in it.product(*tuple(range(-gi, gi+1) for gi in ghosts)):
+ fidx = tuple(np.add(idx, offset))
+ sidx = tuple(np.add(offset, ghosts))
eq -= fvals[fidx]*stencil[sidx]
eqs.append(eq)
@@ -434,7 +470,7 @@ class PolynomialInterpolator(object):
print '\nBuilding matrix...'
M = Ainv.dot(B)
self.M = M
- update_cache(self.cache_file(), self.key, M)
+ update_cache(self.cache_file(), self.key, (P0,S,M))
def interpolate(self, fvals):
@@ -458,7 +494,7 @@ class PolynomialSubgridInterpolator(object):
----------
interpolator: PolynomialInterpolator
Interpolant used to compute weights.
- grid_ratio: tuple of integers
+ grid_ratio: tuple of int
Tuple of integers representing the ratio between the coarse and the fine grid.
dtype: np.dtype
Force to cast dtype for all matrices (interpolator.M may contain rationals).
@@ -467,14 +503,15 @@ class PolynomialSubgridInterpolator(object):
----------
dim: int
Number of dimensions to interpolate (same as interpolator.dim).
- ghosts: int
+ ghosts: tuple of int
Number of required ghosts.
- n: int
+ n: tuple of int
Corresponds to 2*(ghosts+1), the number of required nodes to generate the
polynomial coefficients (same as interpolator.n).
N: int
Total number of input nodes N including ghosts (same as interpolator.N).
- s: int
+ N = n0*n1*...*n[dim-1]
+ s: tuple of int
Corresponds to grid_ratio + 1, number of points of the subgrid in each directions.
Example for a grid ratio=(3,3), we have s=(4,4): O=coarse grid nodes, X=fine grid nodes
@@ -489,10 +526,48 @@ class PolynomialSubgridInterpolator(object):
S: int
Represents the number of fine grid points contained in a coarse grid cell.
S = s0*s1*...*s[dim-1]
+ gr: tuple of int
+ Corresponds to grid_ratio, number of points of the subgrid in each directions, minus one.
+ Example for a grid ratio=(3,3), we have gr=(3,3) and s=(4,4):
+ O=coarse grid nodes, X=fine grid nodes, -=excluded find grid nodes
+
+ Coarse grid: Fine grid:
+
+ ^ O-----O ^ O X X O ^
+ | | | gr0 | X X X - | s0
+ 1 | | | v X X X - |
+ v O-----O O - - O v
+ <-----> <--->
+ 1 gr1
+ GR: int
+ Represents the number of fine grid points contained in a coarse grid cell exluding right most points.
+ GR = gr0*gr1*...*gr[dim-1]
W: np.ndarray
- Pre computed weights to interpolate directly from coarse to fine grid.
+ Pre computed weights to interpolate directly from coarse to fine grid.
+ Let F be the vector of N known coarse grid node values (including required ghosts).
+ Let G be the vector of S unknown fine grid node values.
+
+ N
+ <--------->
+ X X X X X X ^ |g0| ^
+ X X X X X X | |f0| ^ |g1| |
+ X X X X X X | |f1| | |g2| |
+ X X X X X X | |f2| | |g3| |
+ W = X X X X X X | S F= |f3| | N G = W*F = |g4| | S
+ X X X X X X | |f4| | |g5| |
+ X X X X X X | |f5| v |g6| |
+ X X X X X X | |g7| |
+ X X X X X X v |g8| v
+
+ Will contain the same data type as intepolator.M if dtype is not passed,
+ else W will be computed from user given dtype.
+
+ Wr: np.ndarray
+ Reduced W that exludes rightmost output points of ndimensional output vector.
+
+ Pre computed weights to interpolate directly from coarse to inner fine grid.
Let F be the vector of N known coarse grid node values (including required ghosts).
- Let G be the vector of S unknown fine grid node values.
+ Let G be the vector of GR unknown fine inner grid node values (see gr attribute).
N
<--------->
@@ -500,33 +575,33 @@ class PolynomialSubgridInterpolator(object):
X X X X X X | |f0| ^ |g1| |
X X X X X X | |f1| | |g2| |
X X X X X X | |f2| | |g3| |
- W = X X X X X X | S F= |f3| | N G = W*F = |g3| | S
- X X X X X X | |f4| | |g3| |
- X X X X X X | |f5| v |g3| |
- X X X X X X | |g3| |
- X X X X X X v |g4| v
-
- Will contain the same data type as intepolator.M if dtype is not passed,
- else W will be computed from user given dtype.
-
- U: tuple
- Contains unique values contained in W.
- I: np.ndarray of np.int8 or np.int16
- Represents W as integers indexed into U.
+ Wr = X X X X X X | GR F= |f3| | N G = W*F = |g4| | GR
+ X X X X X X | |f4| | |g5| |
+ X X X X X X | |f5| v |g6| |
+ X X X X X X | |g7| |
+ X X X X X X v |g8| v
+
+ Same data type as W.
"""
check_instance(interpolator, PolynomialInterpolator)
check_instance(grid_ratio, tuple, size=interpolator.dim, minval=1)
p = interpolator.p
- P = interpolator.P
n = interpolator.n
+ P = interpolator.P
N = interpolator.N
+ M = interpolator.M
ghosts = interpolator.ghosts
+
+ gr = grid_ratio
+ GR = np.prod(gr, dtype=np.int32)
+ del grid_ratio
- s = tuple(gr+1 for gr in grid_ratio)
- S = np.prod(s, dtype=np.int64)
+ s = tuple(gri+1 for gri in gr)
+ S = np.prod(s, dtype=np.int32)
+
dim = interpolator.dim
- key = ('PolynomialSubgridInterpolator', interpolator.key, grid_ratio, str(dtype))
+ key = ('PolynomialSubgridInterpolator', interpolator.key, gr, str(dtype))
self.p = p
self.P = p
@@ -534,60 +609,48 @@ class PolynomialSubgridInterpolator(object):
self.S = S
self.n = n
self.N = N
+ self.gr = gr
+ self.GR = GR
self.ghosts = interpolator.ghosts
self.dim = dim
self.interpolator = interpolator
self.key = key
cache_file = interpolator.cache_file()
- data = load_data_from_cache(cache_file, key)
- if (data is not None) and _check_matrices(*data):
- (W,U,I) = data
- else:
- X = tuple(np.asarray(tuple(sm.Rational(j,gr) for j in xrange(0,si))) for i,(gr,si) in enumerate(zip(grid_ratio, s)))
- M = interpolator.M
- V = np.vander(X[0], N=p, increasing=True)
- for i in xrange(1, dim):
- Vi = np.vander(X[i], N=p, increasing=True)
- V = np.multiply.outer(V,Vi)
-
- even_axes = range(0,V.ndim,2)
- odd_axes = range(1,V.ndim,2)
- axes = even_axes + odd_axes
-
- V = np.transpose(V, axes=axes).copy()
- assert V.shape[:dim] == s
- assert V.shape[dim:] == (p,)*dim
- V = V.reshape((S,P))
- W = V.dot(M)
- assert W.shape == (S,N)
-
- if interpolator.approximative:
- U = None
- I = None
- else:
- U = list(sorted(set(W.ravel().tolist())))
- U.remove(0)
- U.remove(1)
- U = (0,1,)+tuple(U)
-
- if len(U)<np.iinfo(np.uint8).max:
- itype = np.uint8
- elif len(U)<np.iinfo(np.uint16).max:
- itype = np.uint16
- else:
- raise NotImplementedError(len(U))
-
- I = np.empty(shape=W.shape, dtype=itype)
- for idx in np.ndindex(*W.shape):
- I[idx] = U.index(W[idx])
-
- if (dtype is not None):
- W = W.astype(dtype)
- update_cache(cache_file, key, (W,U,I))
+ try:
+ data = load_data_from_cache(cache_file, key)
+ if (data is not None):
+ W = data
+ if _check_matrix(W):
+ self.W = data
+ return
+ except Exception as e:
+ msg='Failed to load data from cache because:\n{}'.format(e)
+ warnings.warn(msg, HysopCacheWarning)
+
+ X = tuple(np.asarray(tuple(sm.Rational(j,gr) for j in xrange(0,si)))
+ for i,(gr,si) in enumerate(zip(gr, s)))
+ V = np.vander(X[0], N=p[0], increasing=True)
+ for i in xrange(1, dim):
+ Vi = np.vander(X[i], N=p[i], increasing=True)
+ V = np.multiply.outer(V,Vi)
+
+ even_axes = range(0,V.ndim,2)
+ odd_axes = range(1,V.ndim,2)
+ axes = even_axes + odd_axes
+
+ V = np.transpose(V, axes=axes).copy()
+ assert V.shape[:dim] == s
+ assert V.shape[dim:] == p
+ V = V.reshape((S,P))
+ W = V.dot(M)
+ assert W.shape == (S,N)
+
+ if (dtype is not None):
+ W = W.astype(dtype)
+
+ update_cache(cache_file, key, W)
self.W = W
- self.U = U
- self.I = I
def __call__(self, F):
return self.W.dot(F.ravel()).reshape(self.s)
@@ -595,6 +658,13 @@ class PolynomialSubgridInterpolator(object):
def generate_subgrid_restrictor(self):
return PolynomialSubgridRestrictor(subgrid_interpolator=self)
+ @property
+ def Wr(self):
+ assert self.W.shape == (self.S, self.N)
+ view = (slice(0,-1,None),)*self.dim + (slice(None, None, None),)*self.dim
+ Wr = self.W.reshape(self.s+self.n)[view].reshape(self.GR, self.N).copy()
+ return Wr
+
class PolynomialSubgridRestrictor(object):
def __init__(self, subgrid_interpolator):
@@ -677,6 +747,31 @@ if __name__ == '__main__':
print 'Biquintic (FDC2)'
print GI1(F)
print
+
+ F = [[0,1,1,0],
+ [0,1,1,0]]
+ F = np.asarray(F)
+
+ print 'Solving linear/cubic...'
+ PI0 = PolynomialInterpolator(dim=2, deg=(1,3), fd=2, verbose=False)
+ GI0 = PI0.generate_subgrid_interpolator(grid_ratio=grid_ratio, dtype=np.float64)
+ print 'Linear/Cubic (FDC2)'
+ print GI0(F)
+ print
+
+ F = [[0,0],
+ [1,1],
+ [1,1],
+ [0,0]]
+ F = np.asarray(F)
+
+ print 'Solving cubic/linear...'
+ PI0 = PolynomialInterpolator(dim=2, deg=(3,1), fd=2, verbose=False)
+ GI0 = PI0.generate_subgrid_interpolator(grid_ratio=grid_ratio, dtype=np.float64)
+ print 'Cubic/Linear (FDC2)'
+ print GI0(F)
+ print
+
F = [[0,0,0,0,0,0],
[0,0,0,0,0,0],
@@ -714,6 +809,20 @@ if __name__ == '__main__':
print GI3(F)
print
+ print 'Solving septic2/nonic2 ...'
+ PI4 = PolynomialInterpolator(dim=2, deg=(7,9), fd=2, verbose=False, approximative=(not has_flint))
+ GI4 = PI4.generate_subgrid_interpolator(grid_ratio=grid_ratio, dtype=np.float64)
+ print 'septic/nonic (FDC2)'
+ print GI4(F)
+ print
+
+ print 'Solving septic2/quintic4 ...'
+ PI5 = PolynomialInterpolator(dim=2, deg=(7,5), fd=(2,4), verbose=False, approximative=(not has_flint))
+ GI5 = PI5.generate_subgrid_interpolator(grid_ratio=grid_ratio, dtype=np.float64)
+ print 'septic/nonic (FDC2/FDC4)'
+ print GI5(F)
+ print
+
# 3D test
grid_ratio = (2,2,2)
print 'Solving trilinear...'
diff --git a/hysop/operator/base/spatial_filtering.py b/hysop/operator/base/spatial_filtering.py
index 401a02e1c51291318d120876447ed5e5affff740..a73b8e975113cbcd89644d4360d11e87ad72f855 100644
--- a/hysop/operator/base/spatial_filtering.py
+++ b/hysop/operator/base/spatial_filtering.py
@@ -2,6 +2,7 @@
@file advection.py
RestrictionFilter operator generator.
"""
+import numpy as np
from hysop.constants import Implementation
from hysop.methods import Remesh
from hysop.numerics.remesh.remesh import RemeshKernel
@@ -39,11 +40,17 @@ class SpatialFilterBase(object):
input_fields={input_field: input_topo},
output_fields={output_field: output_topo},
**kwds)
- assert (input_field.dim == output_field.dim)
- self.Fin = input_field
- self.Fout = output_field
- self.dim = input_field.dim
- self.dtype = find_common_dtype(input_field.dtype, output_field.dtype)
+
+ Fin = input_field
+ Fout = output_field
+ assert (Fin.dim == Fout.dim)
+ assert (Fin.lboundaries == Fout.lboundaries).all()
+ assert (Fin.rboundaries == Fout.rboundaries).all()
+ assert (Fin.periodicity == Fout.periodicity).all()
+ self.Fin = Fin
+ self.Fout = Fout
+ self.dim = Fin.dim
+ self.dtype = find_common_dtype(Fin.dtype, Fout.dtype)
self.iratio = None # will be set in get_field_requirements
self.grid_ratio = None # will be set in discretize
@@ -58,6 +65,7 @@ class SpatialFilterBase(object):
self.dFin = dFin
self.dFout = dFout
self.grid_ratio = grid_ratio
+ self.iter_shape = self.dFin.compute_resolution + 1 - self.dFin.periodicity
@classmethod
def supports_multiple_field_topologies(cls):
@@ -381,7 +389,7 @@ class PolynomialInterpolationFilterBase(PolynomialInterpolationMethod, Interpola
@debug
def get_field_requirements(self):
reqs = super(PolynomialInterpolationFilterBase, self).get_field_requirements()
- required_input_ghosts = (self.polynomial_interpolator.ghosts+1,)*self.dim
+ required_input_ghosts = np.add(self.polynomial_interpolator.ghosts, self.Fin.periodicity)
Fin_topo, Fin_requirements = reqs.get_input_requirement(self.Fin)
Fin_requirements.min_ghosts = required_input_ghosts
self.required_input_ghosts = required_input_ghosts
@@ -392,8 +400,9 @@ class PolynomialInterpolationFilterBase(PolynomialInterpolationMethod, Interpola
return
super(PolynomialInterpolationFilterBase, self).discretize()
dFin, dFout = self.dFin, self.dFout
+ ghosts = self.dFin.topology_state.transposed(self.required_input_ghosts)
psi = self.polynomial_interpolator.generate_subgrid_interpolator(grid_ratio=self.grid_ratio)
self.subgrid_interpolator = psi
- self.fin = dFin.sdata[dFin.local_slices(ghosts=self.required_input_ghosts)].handle
+ self.fin = dFin.sdata[dFin.local_slices(ghosts=ghosts)].handle
self.fout = dFout.sdata[dFout.compute_slices].handle
diff --git a/hysop/tools/warning.py b/hysop/tools/warning.py
index 18da7991bf9052b22bed8b5c2ebf7bb06329e547..cb040bf5b18d67b0b8e576daff91d4ad4e0cd136 100644
--- a/hysop/tools/warning.py
+++ b/hysop/tools/warning.py
@@ -17,6 +17,12 @@ class HysopPerformanceWarning(HysopWarning):
"""
pass
+class HysopCacheWarning(HysopWarning):
+ """
+ Custom warning class for hysop caching.
+ """
+ pass
+
def configure_hysop_warnings(action):
"""
Configure hysop warnings.