From 074350c93adbf3070cfc329ebf344b70fdee3f74 Mon Sep 17 00:00:00 2001
From: Jean-Baptiste Keck <Jean-Baptiste.Keck@imag.fr>
Date: Fri, 12 Apr 2019 11:05:29 +0200
Subject: [PATCH] working polynomial filters
---
examples/multiresolution/scalar_advection.py | 15 ++-
.../opencl/operator/spatial_filtering.py | 52 +++++++-
.../host/python/operator/spatial_filtering.py | 38 +++++-
hysop/numerics/interpolation/polynomial.py | 111 +++++++++++++-----
hysop/operator/base/spatial_filtering.py | 35 +++++-
hysop/operator/spatial_filtering.py | 29 +++--
hysop/symbolic/array.py | 3 +-
7 files changed, 229 insertions(+), 54 deletions(-)
diff --git a/examples/multiresolution/scalar_advection.py b/examples/multiresolution/scalar_advection.py
index 885d756dd..6b68a8a7e 100644
--- a/examples/multiresolution/scalar_advection.py
+++ b/examples/multiresolution/scalar_advection.py
@@ -30,6 +30,7 @@ def compute(args):
npts = args.npts # coarse resolution
snpts = args.snpts # fine resolution
fnpts = tuple(3*_ for _ in snpts) # finest resolution
+ cnpts = tuple(_//2 for _ in npts) # coarsest resolution
lboundaries = (BoxBoundaryCondition.PERIODIC,)*dim
rboundaries = (BoxBoundaryCondition.PERIODIC,)*dim
box = Box(origin=args.box_origin, length=args.box_length, dim=dim,
@@ -118,9 +119,9 @@ def compute(args):
filtering_method=args.interpolation_filter,
implementation=impl,
**extra_op_kwds)
- #> Restriction filter
- restriction_filter = SpatialFilter(input_variables={scalar: snpts},
- output_variables={scalar: npts},
+ #> Restriction filter
+ restriction_filter = SpatialFilter(input_variables={scalar: npts},
+ output_variables={scalar: cnpts},
filtering_method=args.restriction_filter,
implementation=impl,
**extra_op_kwds)
@@ -141,9 +142,15 @@ def compute(args):
io_params=io_params,
variables={scalar: npts},
**extra_op_kwds)
+ io_params = IOParams(filename='coarsest', frequency=args.dump_freq)
+ df3 = HDF_Writer(name='S_coarsest',
+ io_params=io_params,
+ variables={scalar: cnpts},
+ **extra_op_kwds)
# Add a writer of input field at given frequency.
- problem.insert(interpolation_filter, df0, df1, df2)
+ problem.insert(interpolation_filter, restriction_filter,
+ df0, df1, df2, df3)
problem.build(args)
# If a visu_rank was provided, and show_graph was set,
diff --git a/hysop/backend/device/opencl/operator/spatial_filtering.py b/hysop/backend/device/opencl/operator/spatial_filtering.py
index 0b4271c79..1501c1b29 100644
--- a/hysop/backend/device/opencl/operator/spatial_filtering.py
+++ b/hysop/backend/device/opencl/operator/spatial_filtering.py
@@ -11,7 +11,7 @@ from hysop.fields.continuous_field import Field
from hysop.parameters.parameter import Parameter
from hysop.topology.cartesian_descriptor import CartesianTopologyDescriptors
from hysop.operator.base.spatial_filtering import RemeshRestrictionFilterBase, SpectralRestrictionFilterBase, \
- SubgridRestrictionFilterBase, PolynomialInterpolationFilterBase
+ SubgridRestrictionFilterBase, PolynomialInterpolationFilterBase, PolynomialRestrictionFilterBase
from hysop.backend.device.opencl.opencl_symbolic import OpenClSymbolic
from hysop.backend.device.opencl.opencl_copy_kernel_launchers import OpenClCopyBufferRectLauncher
from hysop.backend.device.opencl.opencl_kernel_launcher import OpenClKernelListLauncher
@@ -57,7 +57,7 @@ class OpenClPolynomialInterpolationFilter(PolynomialInterpolationFilterBase, Ope
exchange_ghosts = self.dFout.exchange_ghosts(build_launcher=True)
- kl = OpenClKernelListLauncher(name='lowpass_filter')
+ kl = OpenClKernelListLauncher(name=kname)
kl += interpolate_grid_kernel
kl += exchange_ghosts
@@ -69,6 +69,54 @@ class OpenClPolynomialInterpolationFilter(PolynomialInterpolationFilterBase, Ope
evt = self.execute_kernels()
+class OpenClPolynomialRestrictionFilter(PolynomialRestrictionFilterBase, OpenClOperator):
+
+ @debug
+ def discretize(self):
+ if self.discretized:
+ return
+ super(OpenClPolynomialRestrictionFilter, self).discretize()
+ dFin = self.dFin
+ dFout = self.dFout
+ gr = self.grid_ratio
+ dim = dFin.dim
+ assert dFin.is_scalar
+ assert dFout.is_scalar
+ assert self.subgrid_restrictor.gr == gr
+
+ ekg = self.elementwise_kernel_generator
+ Rr = self.subgrid_restrictor.Rr / self.subgrid_restrictor.GR
+ ghosts = np.asarray(self.subgrid_restrictor.ghosts)
+
+ I = np.asarray(local_indices_symbols[:dim][::-1])
+ fin, fout = ekg.dfields_to_ndbuffers(dFin, dFout)
+
+ def gen_inputs(*idx):
+ return fin(gr*I+idx-ghosts)
+ input_values = np.asarray(tuple(map(gen_inputs, np.ndindex(*Rr.shape)))).reshape(Rr.shape)
+ output_value = (Rr*input_values).sum()
+
+ e = Assignment(fout(I), output_value)
+ exprs = (e,)
+
+ kname='restrict_grid_{}'.format(self.polynomial_interpolation_method).lower()
+ restriction_grid_kernel, _ = ekg.elementwise_kernel(kname,
+ *exprs, compute_resolution=self.iter_shape, debug=False)
+
+ exchange_ghosts = self.dFout.exchange_ghosts(build_launcher=True)
+
+ kl = OpenClKernelListLauncher(name=kname)
+ kl += restriction_grid_kernel
+ kl += exchange_ghosts
+
+ self.execute_kernels = functools.partial(kl, queue=self.cl_env.default_queue)
+
+ @op_apply
+ def apply(self, **kwds):
+ super(OpenClPolynomialRestrictionFilter, self).apply(**kwds)
+ evt = self.execute_kernels()
+
+
class OpenClSubgridRestrictionFilter(SubgridRestrictionFilterBase, OpenClSymbolic):
"""
OpenCL implementation for lowpass spatial filtering: small grid -> coarse grid
diff --git a/hysop/backend/host/python/operator/spatial_filtering.py b/hysop/backend/host/python/operator/spatial_filtering.py
index a82688827..811384d31 100644
--- a/hysop/backend/host/python/operator/spatial_filtering.py
+++ b/hysop/backend/host/python/operator/spatial_filtering.py
@@ -7,8 +7,12 @@ from hysop.core.graph.graph import op_apply
from hysop.fields.continuous_field import Field
from hysop.parameters.parameter import Parameter
from hysop.topology.cartesian_descriptor import CartesianTopologyDescriptors
-from hysop.operator.base.spatial_filtering import RemeshRestrictionFilterBase, SpectralRestrictionFilterBase, \
- SubgridRestrictionFilterBase, PolynomialInterpolationFilterBase
+from hysop.operator.base.spatial_filtering import (
+ PolynomialInterpolationFilterBase,
+ RemeshRestrictionFilterBase,
+ SpectralRestrictionFilterBase,
+ SubgridRestrictionFilterBase,
+ PolynomialRestrictionFilterBase)
class PythonPolynomialInterpolationFilter(PolynomialInterpolationFilterBase, HostOperator):
@@ -24,20 +28,42 @@ class PythonPolynomialInterpolationFilter(PolynomialInterpolationFilterBase, Hos
super(PythonPolynomialInterpolationFilter, self).apply(**kwds)
fin = self.fin
fout = self.fout
- ishape = self.dFin.compute_resolution
- oshape = self.dFout.compute_resolution
periodicity = self.dFin.periodicity
- gr = self.grid_ratio
gr, n = self.subgrid_interpolator.gr, self.subgrid_interpolator.n
Wr = self.Wr
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))
+ 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()
+class PythonPolynomialRestrictionFilter(PolynomialRestrictionFilterBase, HostOperator):
+ def discretize(self, **kwds):
+ if self.discretized:
+ return
+ super(PythonPolynomialRestrictionFilter, self).discretize(**kwds)
+ SR = self.subgrid_restrictor
+ self.Rr = SR.Rr.astype(self.dtype) / SR.GR
+ assert (self.Rr.shape == tuple(2*gi+1 for gi in SR.ghosts)), Rr.shape
+
+ @op_apply
+ def apply(self, **kwds):
+ """Apply analytic formula."""
+ super(PythonPolynomialRestrictionFilter, self).apply(**kwds)
+ fin = self.fin
+ fout = self.fout
+ gr = self.subgrid_restrictor.gr
+ Rr = self.Rr
+ rshape = Rr.shape
+
+ for idx in np.ndindex(*self.iter_shape):
+ islc = tuple(slice(j*gr[i], j*gr[i]+rshape[i], 1) for i,j in enumerate(idx))
+ fout[idx] = (Rr*fin[islc]).sum()
+ self.dFout.exchange_ghosts()
+
class PythonRemeshRestrictionFilter(RemeshRestrictionFilterBase, HostOperator):
"""
Python implementation for lowpass spatial filtering: small grid -> coarse grid
diff --git a/hysop/numerics/interpolation/polynomial.py b/hysop/numerics/interpolation/polynomial.py
index b1a27c2aa..70e83e44b 100644
--- a/hysop/numerics/interpolation/polynomial.py
+++ b/hysop/numerics/interpolation/polynomial.py
@@ -116,14 +116,16 @@ class PolynomialInterpolator(object):
Polynomial degree (1=linear, 3=cubic, 5=quintic, 7=septic, ...)
Degree should be odd: deg=2k+1
fd: tuple of ints
- Order of centered finite differences stencils used to compute derivatives for each direction.
+ 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=p0*p1*...*p(dim-1).
P: int
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.
+ 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: tuple of ints
Return the number of ghosts required by the interpolator on each axis.
@@ -166,7 +168,8 @@ class PolynomialInterpolator(object):
See Also
--------
- :class:`PolynomialSubgridInterpolator`: Precompute weights for fixed subgrid interpolation.
+ :class:`PolynomialSubgridInterpolator`: Precompute weights for fixed subgrid
+ interpolation.
"""
assert dim>0, 'dim<=0'
@@ -365,7 +368,8 @@ class PolynomialInterpolator(object):
if verbose:
print ' {}'.format(eq)
- # Build system such that A*c = B*f where c are the polynomial coefficients and f the node values
+ # Build system such that A*c = B*f where c are the polynomial coefficients and
+ # f the node values
dtype = (np.float64 if approximative else object)
A = np.empty((self.P,self.P), dtype=dtype)
B = np.empty((self.P,self.N), dtype=dtype)
@@ -408,13 +412,15 @@ class PolynomialInterpolator(object):
return sm.utilities.lambdify(self.xvars, P0)
def generate_subgrid_interpolator(self, grid_ratio, dtype=None):
- return PolynomialSubgridInterpolator(interpolator=self, grid_ratio=grid_ratio, dtype=dtype)
+ return PolynomialSubgridInterpolator(interpolator=self,
+ grid_ratio=grid_ratio, dtype=dtype)
class PolynomialSubgridInterpolator(object):
def __init__(self, interpolator, grid_ratio, dtype=None):
"""
- Create a PolynomialSubgridInterpolator from a PolynomialInterpolator and a number of subrid points.
+ Create a PolynomialSubgridInterpolator from a PolynomialInterpolator and a number of
+ subrid points.
Parameters
----------
@@ -439,7 +445,8 @@ class PolynomialSubgridInterpolator(object):
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
+ Example for a grid ratio=(3,3), we have s=(4,4):
+ O=coarse grid nodes, X=fine grid nodes
Coarse grid: Fine grid:
@@ -453,7 +460,8 @@ class PolynomialSubgridInterpolator(object):
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.
+ 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
@@ -466,11 +474,13 @@ class PolynomialSubgridInterpolator(object):
<-----> <--->
1 gr1
GR: int
- Represents the number of fine grid points contained in a coarse grid cell exluding right most points.
+ 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.
- Let F be the vector of N known coarse grid node values (including required ghosts).
+ 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
@@ -604,39 +614,63 @@ class PolynomialSubgridRestrictor(object):
Attributes
----------
- g: int
+ g: tuple of int
Corresponds to (n+1)*s
G: int
Corresponds to g[0]*g[1]*...g[dim-1]
+ R: np.ndarray
+ Restrictor weights.
+ origin: tuple of int
+ Origin of the generated stencil.
+ Rr: np.ndarray
+ Restrictor weights excluding leftmost and rightmost points.
+ ghosts:
+ Corresponds to origin - 1, which is also Rr origin.
"""
check_instance(subgrid_interpolator, PolynomialSubgridInterpolator)
- s = subgrid_interpolator.s
- n = subgrid_interpolator.n
- gr = tuple(si-1 for si in s)
- g = tuple(n*gri+1 for gri in gr)
+ dim = subgrid_interpolator.dim
+ n = subgrid_interpolator.n
+ s = subgrid_interpolator.s
+ gr = subgrid_interpolator.gr
+ GR = subgrid_interpolator.GR
+ W = subgrid_interpolator.W
+ g = tuple(ni*gri+1 for (ni,gri) in zip(n,gr))
+ G = np.prod(g, dtype=np.int64)
assert all(gi%2==1 for gi in g)
origin = tuple(gi/2 for gi in g)
- G = np.prod(g, dtype=np.int64)
-
gvals, gvars = tensor_symbol('g',g,origin)
I = 0
for idx in np.ndindex(*gr):
- mask = tuple(slice(i, max(2,i+n*gri), gri) for (i,gri) in zip(idx, gr))
+ mask = tuple(slice(i, max(2,i+ni*gri), gri) for (i,ni,gri) in zip(idx, n, gr))
target = tuple(gri-i for (i,gri) in zip(idx,gr))
F = gvals[mask]
- Ii = subgrid_interpolator.W.dot(F.ravel()).reshape(s)[target]
+ Ii = W.dot(F.ravel()).reshape(s)[target]
I += Ii
R = np.ndarray(shape=g, dtype=object)
for idx in np.ndindex(*g):
R[idx] = I.coeff(gvals[idx])
- print R.astype(np.float32)
- print R.sum()
+
+ view = (slice(1,-1,None),)*dim
+ Rr = R[view]
+ ghosts = tuple(oi-1 for oi in origin)
+
+ self.g = g
+ self.G = G
+ self.R = R
+ self.Rr = Rr
+ self.origin = origin
+ self.ghosts = ghosts
+ self.n = n
+ self.s = s
+ self.GR = GR
+ self.gr = gr
+ self.subgrid_interpolator = subgrid_interpolator
if __name__ == '__main__':
np.set_printoptions(precision=4, linewidth=1e8, threshold=1e8,
formatter={'float': lambda x: "{0:+0.3f}".format(x)})
-
+
# 2D tests
grid_ratio = (2,2)
F = [[1,1],
@@ -722,28 +756,32 @@ if __name__ == '__main__':
print
print 'Solving biseptic2...'
- PI2 = PolynomialInterpolator(dim=2, deg=7, fd=2, verbose=False, approximative=(not has_flint))
+ PI2 = PolynomialInterpolator(dim=2, deg=7, fd=2, verbose=False,
+ approximative=(not has_flint))
GI2 = PI2.generate_subgrid_interpolator(grid_ratio=grid_ratio, dtype=np.float64)
print 'Biseptic (FDC2)'
print GI2(F)
print
print 'Solving binonic2...'
- PI3 = PolynomialInterpolator(dim=2, deg=9, fd=2, verbose=False, approximative=(not has_flint))
+ PI3 = PolynomialInterpolator(dim=2, deg=9, fd=2, verbose=False,
+ approximative=(not has_flint))
GI3 = PI3.generate_subgrid_interpolator(grid_ratio=grid_ratio, dtype=np.float64)
print 'Binonic (FDC2)'
print GI3(F)
print
print 'Solving septic2/nonic2 ...'
- PI4 = PolynomialInterpolator(dim=2, deg=(7,9), fd=2, verbose=False, approximative=(not has_flint))
+ 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))
+ 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)
@@ -762,8 +800,25 @@ if __name__ == '__main__':
print
print 'Solving triquintic2...'
- PI0 = PolynomialInterpolator(dim=3, deg=5, fd=2, verbose=False, approximative=(not has_flint))
+ PI0 = PolynomialInterpolator(dim=3, deg=5, fd=2, verbose=False,
+ approximative=(not has_flint))
GI0 = PI0.generate_subgrid_interpolator(grid_ratio=grid_ratio, dtype=np.float64)
print
-
+ # Delaurier-Dubuc interpolating wavelets
+ from matplotlib import pyplot as plt
+ grid_ratio = (32,)
+ fig, axes = plt.subplots(ncols=2, nrows=2)
+ #for k,(deg, sdeg) in enumerate(zip((3,5,7,9),('cubic','quintic','septic','nonic'))):
+ for k,(deg, sdeg) in enumerate(zip((3,5,),('cubic','quintic'))):
+ ax = axes[k/2, k%2]
+ for fd in (2,4,6):
+ PI = PolynomialInterpolator(dim=1, deg=deg, fd=fd, verbose=False)
+ SI = PI.generate_subgrid_interpolator(grid_ratio=grid_ratio)
+ SR = SI.generate_subgrid_restrictor()
+ X = np.linspace(-SR.n[0]/2, +SR.n[0]/2, SR.g[0])
+ ax.plot(X, SR.R, label='{}_fdc{}'.format(sdeg, fd))
+ ax.legend()
+ ax.plot(X, np.zeros(SR.g[0]), '--')
+ plt.show()
+
diff --git a/hysop/operator/base/spatial_filtering.py b/hysop/operator/base/spatial_filtering.py
index 793fb86e4..30db42947 100644
--- a/hysop/operator/base/spatial_filtering.py
+++ b/hysop/operator/base/spatial_filtering.py
@@ -65,7 +65,6 @@ 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):
@@ -401,8 +400,40 @@ class PolynomialInterpolationFilterBase(PolynomialInterpolationMethod, Interpola
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)
+ psi = self.polynomial_interpolator.generate_subgrid_interpolator(
+ grid_ratio=self.grid_ratio)
self.subgrid_interpolator = psi
self.fin = dFin.sdata[dFin.local_slices(ghosts=ghosts)].handle
self.fout = dFout.sdata[dFout.compute_slices].handle
+ self.iter_shape = self.dFin.compute_resolution + 1 - self.dFin.periodicity
+
+
+class PolynomialRestrictionFilterBase(PolynomialInterpolationMethod, RestrictionFilterBase):
+ """
+ Base implementation for polynomial interpolation.
+ """
+ @debug
+ def get_field_requirements(self):
+ reqs = super(PolynomialRestrictionFilterBase, self).get_field_requirements()
+ iratio = self.iratio
+ pghosts = self.polynomial_interpolator.ghosts
+ ghosts = np.add(np.multiply(iratio, np.add(pghosts,1)), -1)
+ Fin_topo, Fin_requirements = reqs.get_input_requirement(self.Fin)
+ Fin_requirements.min_ghosts = ghosts
+ self.required_input_ghosts = ghosts
+ return reqs
+
+ def discretize(self):
+ if self.discretized:
+ return
+ super(PolynomialRestrictionFilterBase, self).discretize()
+ dFin, dFout = self.dFin, self.dFout
+ ghosts = self.dFin.topology_state.transposed(self.required_input_ghosts)
+ psr = self.polynomial_interpolator.generate_subgrid_interpolator(
+ grid_ratio=self.grid_ratio).generate_subgrid_restrictor()
+ assert all(psr.ghosts == ghosts)
+ self.subgrid_restrictor = psr
+ self.fin = dFin.sdata[dFin.local_slices(ghosts=ghosts)].handle
+ self.fout = dFout.sdata[dFout.compute_slices].handle
+ self.iter_shape = self.dFout.compute_resolution
diff --git a/hysop/operator/spatial_filtering.py b/hysop/operator/spatial_filtering.py
index 38d2f4e58..f485aec69 100644
--- a/hysop/operator/spatial_filtering.py
+++ b/hysop/operator/spatial_filtering.py
@@ -77,20 +77,26 @@ class RestrictionFilterFrontend(SpatialFilterFrontend):
@classmethod
def all_implementations(cls):
from hysop.backend.host.python.operator.spatial_filtering import \
- PythonRemeshRestrictionFilter, PythonSpectralRestrictionFilter, PythonSubgridRestrictionFilter
+ PythonRemeshRestrictionFilter, PythonSpectralRestrictionFilter, \
+ PythonSubgridRestrictionFilter, PythonPolynomialRestrictionFilter
from hysop.backend.device.opencl.operator.spatial_filtering import \
- OpenClSpectralRestrictionFilter, OpenClSubgridRestrictionFilter
+ OpenClSpectralRestrictionFilter, OpenClSubgridRestrictionFilter, \
+ OpenClPolynomialRestrictionFilter
ai = {
+ FilteringMethod.SUBGRID: {
+ Implementation.PYTHON: PythonSubgridRestrictionFilter,
+ Implementation.OPENCL: OpenClSubgridRestrictionFilter,
+ },
+ FilteringMethod.POLYNOMIAL: {
+ Implementation.PYTHON: PythonPolynomialRestrictionFilter,
+ Implementation.OPENCL: OpenClPolynomialRestrictionFilter,
+ },
FilteringMethod.SPECTRAL: {
Implementation.PYTHON: PythonSpectralRestrictionFilter,
- Implementation.OPENCL: OpenClSpectralRestrictionFilter
+ Implementation.OPENCL: OpenClSpectralRestrictionFilter,
},
FilteringMethod.REMESH: {
- Implementation.PYTHON: PythonRemeshRestrictionFilter
- },
- FilteringMethod.SUBGRID: {
- Implementation.PYTHON: PythonSubgridRestrictionFilter,
- Implementation.OPENCL: OpenClSubgridRestrictionFilter
+ Implementation.PYTHON: PythonRemeshRestrictionFilter,
},
}
return ai
@@ -98,9 +104,10 @@ class RestrictionFilterFrontend(SpatialFilterFrontend):
@classmethod
def all_default_implementations(cls):
adi = {
- FilteringMethod.REMESH: Implementation.PYTHON,
- FilteringMethod.SPECTRAL: Implementation.PYTHON,
- FilteringMethod.SUBGRID: Implementation.PYTHON
+ FilteringMethod.SUBGRID: Implementation.PYTHON,
+ FilteringMethod.POLYNOMIAL: Implementation.PYTHON,
+ FilteringMethod.SPECTRAL: Implementation.PYTHON,
+ FilteringMethod.REMESH: Implementation.PYTHON,
}
return adi
diff --git a/hysop/symbolic/array.py b/hysop/symbolic/array.py
index 4c28c0765..a1ae050f0 100644
--- a/hysop/symbolic/array.py
+++ b/hysop/symbolic/array.py
@@ -257,7 +257,8 @@ class SymbolicNdBuffer(SymbolicBuffer):
def __call__(self, *idx):
if len(idx)==1 and isinstance(idx[0], npw.ndarray):
- idx = tuple(idx[0].tolist())
+ assert idx[0].size==self._dim, idx[0].shape
+ idx = tuple(idx[0].ravel().tolist())
assert len(idx) == self._dim, idx
offset = npw.dot(self._symbolic_strides, npw.add(idx, self._symbolic_ghosts))
return self.__getitem__(key=offset)
--
GitLab