From 3fa0e1ba9826cd89909bb757bf27cd9021c33028 Mon Sep 17 00:00:00 2001
From: Jean-Baptiste Keck <Jean-Baptiste.Keck@imag.fr>
Date: Wed, 10 Apr 2019 17:06:25 +0200
Subject: [PATCH] opencl interpolation
---
.../backend/device/opencl/opencl_symbolic.py | 140 +++++-------------
.../opencl/operator/spatial_filtering.py | 28 +++-
hysop/numerics/interpolation/polynomial.py | 6 +-
hysop/symbolic/array.py | 18 ++-
hysop/tools/sympy_utils.py | 36 +++--
5 files changed, 100 insertions(+), 128 deletions(-)
diff --git a/hysop/backend/device/opencl/opencl_symbolic.py b/hysop/backend/device/opencl/opencl_symbolic.py
index f0939dca0..d0d414e62 100644
--- a/hysop/backend/device/opencl/opencl_symbolic.py
+++ b/hysop/backend/device/opencl/opencl_symbolic.py
@@ -9,9 +9,10 @@ expressions.
used to provide a common interface to all discrete operators working with the
opencl backend and using kernels generated on the fly from symbolic expressions.
"""
+import numpy as np
from abc import ABCMeta
from hysop.tools.decorators import debug
-from hysop.tools.types import check_instance, first_not_None, InstanceOf
+from hysop.tools.types import check_instance, first_not_None, InstanceOf, to_tuple
from hysop.tools.numpywrappers import npw
from hysop.fields.continuous_field import ScalarField, TensorField, Field
from hysop.parameters.tensor_parameter import Parameter, TensorParameter
@@ -80,124 +81,65 @@ class OpenClSymbolic(OpenClOperator):
self.extra_kwds[name] = extra_kwds
@classmethod
- def symbolic_ndbuffers(cls, *names, **kwds):
- from hysop.symbolic.array import OpenClSymbolicNdBuffer
- buffers = ()
- if ('count' in kwds):
- count = kwds.pop('count')
+ def __symbolic_variables(cls, *names, **kwds):
+ scls = kwds.pop('scls')
+
+ arrays = ()
+ shape = to_tuple(kwds.get('shape', kwds.get('count', ())))
+ if shape:
assert len(names)==1
- snames = tuple(names[0]+subscript(i) for i in xrange(count))
- varnames = tuple('{}{}'.format(names[0], i) for i in xrange(count))
- latexnames = tuple('{}_{{{}}}'.format(names[0], i) for i in xrange(count))
+ snames = tuple(names[0]+subscripts(idx,',',disable_unicode=False) for idx in np.ndindex(*shape))
+ var_names = tuple(names[0]+subscripts(idx,'_',disable_unicode=True) for idx in np.ndindex(*shape))
else:
- snames = names
- varnames = names
- latexnames = names
- for (name, pname, var_name, latex_name) in zip(varnames, snames, varnames, latexnames):
+ snames = names
+ var_names = names
+ shape = len(names)
+ for (name, pname, var_name) in zip(var_names, snames, var_names):
assert ',' not in name
check_instance(name, str)
- buf = OpenClSymbolicNdBuffer(name=name,
- pretty_name=pname,
- var_name=var_name,
- latex_name=latex_name,
- memory_object=None, **kwds)
- buffers += (buf,)
- return buffers
+ arr = scls(name=name, pretty_name=pname, var_name=var_name, **kwds)
+ arrays += (arr,)
+ return np.asarray(arrays).reshape(shape)
+
+ @classmethod
+ def symbolic_ndbuffers(cls, *names, **kwds):
+ from hysop.symbolic.array import OpenClSymbolicNdBuffer
+ assert 'memory_object' not in kwds
+ assert 'scls' not in kwds
+ kwds['memory_object'] = None
+ kwds['scls'] = OpenClSymbolicNdBuffer
+ return cls.__symbolic_variables(*names, **kwds)
@classmethod
def symbolic_buffers(cls, *names, **kwds):
from hysop.symbolic.array import OpenClSymbolicBuffer
- buffers = ()
- if ('count' in kwds):
- count = kwds.pop('count')
- assert len(names)==1
- snames = tuple(names[0]+subscript(i) for i in xrange(count))
- varnames = tuple('{}{}'.format(names[0], i) for i in xrange(count))
- latexnames = tuple('{}_{{{}}}'.format(names[0], i) for i in xrange(count))
- else:
- snames = names
- varnames = names
- latexnames = names
- for (name, pname, var_name, latex_name) in zip(varnames, snames, varnames, latexnames):
- assert ',' not in name
- check_instance(name, str)
- buf = OpenClSymbolicBuffer(name=name,
- pretty_name=pname,
- var_name=var_name,
- latex_name=latex_name,
- memory_object=None, **kwds)
- buffers += (buf,)
- return buffers
+ assert 'memory_object' not in kwds
+ assert 'scls' not in kwds
+ kwds['memory_object'] = None
+ kwds['scls'] = OpenClSymbolicBuffer
+ return cls.__symbolic_variables(*names, **kwds)
@classmethod
def symbolic_arrays(cls, *names, **kwds):
from hysop.symbolic.array import OpenClSymbolicArray
- arrays = ()
- if ('count' in kwds):
- count = kwds.pop('count')
- assert len(names)==1
- snames = tuple(names[0]+subscript(i) for i in xrange(count))
- var_names = tuple('{}{}'.format(names[0], i) for i in xrange(count))
- latex_names = tuple('{}_{{{}}}'.format(names[0], i) for i in xrange(count))
- else:
- snames = names
- var_names = names
- latex_names = names
- for (name, pname, var_name, latex_name) in zip(var_names, snames, var_names, latex_names):
- assert ',' not in name
- check_instance(name, str)
- arr = OpenClSymbolicArray(name=name, pretty_name=pname,
- var_name=var_name, latex_name=latex_name,
- memory_object=None, **kwds)
- arrays += (arr,)
- return arrays
+ assert 'memory_object' not in kwds
+ assert 'scls' not in kwds
+ kwds['scls'] = OpenClSymbolicArray
+ return cls.__symbolic_variables(*names, **kwds)
@classmethod
def symbolic_tmp_scalars(cls, *names, **kwds):
from hysop.symbolic.tmp import TmpScalar
- scalars = ()
- if ('count' in kwds):
- count = kwds.pop('count')
- assert len(names)==1
- snames = tuple(names[0]+subscript(i) for i in xrange(count))
- var_names = tuple('{}{}'.format(names[0], i) for i in xrange(count))
- latex_names = tuple('{}_{{{}}}'.format(names[0], i) for i in xrange(count))
- else:
- snames = names
- var_names = names
- latex_names = names
- for (name, pname, var_name, latex_name) in zip(var_names, snames, var_names, latex_names):
- assert ',' not in name
- check_instance(name, str)
- tmp = TmpScalar(name=name, pretty_name=pname,
- var_name=var_name, latex_name=latex_name,
- **kwds)
- scalars += (tmp,)
- return scalars
+ assert 'scls' not in kwds
+ kwds['scls'] = TmpScalar
+ return cls.__symbolic_variables(*names, **kwds)
@classmethod
def symbolic_constants(cls, *names, **kwds):
- # constants that may be specified up to the setup step
from hysop.symbolic.constant import SymbolicConstant
- constants = ()
- if ('count' in kwds):
- count = kwds.pop('count')
- assert len(names)==1
- snames = tuple(names[0]+subscript(i) for i in xrange(count))
- var_names = tuple('{}{}'.format(names[0], i) for i in xrange(count))
- latex_names = tuple('{}_{{{}}}'.format(names[0], i) for i in xrange(count))
- else:
- snames = names
- var_names = names
- latex_names = names
- for (name, pname, var_name, latex_name) in zip(var_names, snames, var_names, latex_names):
- assert ',' not in name
- check_instance(name, str)
- tmp = SymbolicConstant(name=name, pretty_name=pname,
- var_name=var_name, latex_name=latex_name,
- **kwds)
- constants += (tmp,)
- return constants
+ assert 'scls' not in kwds
+ kwds['scls'] = SymbolicConstant
+ return cls.__symbolic_variables(*names, **kwds)
@debug
def handle_method(self, method):
diff --git a/hysop/backend/device/opencl/operator/spatial_filtering.py b/hysop/backend/device/opencl/operator/spatial_filtering.py
index 520b9f298..c01c11a2b 100644
--- a/hysop/backend/device/opencl/operator/spatial_filtering.py
+++ b/hysop/backend/device/opencl/operator/spatial_filtering.py
@@ -28,6 +28,7 @@ class OpenClPolynomialInterpolationFilter(PolynomialInterpolationFilterBase, Ope
super(OpenClPolynomialInterpolationFilter, self).discretize()
dFin = self.dFin
dFout = self.dFout
+ gr = self.grid_ratio
dim = dFin.dim
assert dFin.is_scalar
assert dFout.is_scalar
@@ -35,24 +36,37 @@ class OpenClPolynomialInterpolationFilter(PolynomialInterpolationFilterBase, Ope
ekg = self.elementwise_kernel_generator
W = self.subgrid_interpolator.W
s = self.subgrid_interpolator.s
- gr = self.grid_ratio
+ n = (self.subgrid_interpolator.n,)*dim
+ ghosts = np.asarray((self.subgrid_interpolator.ghosts,)*dim)
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)
exprs = ()
- for idx in np.ndindex(*s):
- e = Assignment(fin(I+idx), fout(I+idx))
+ for idx in np.ndindex(*n):
+ e = Assignment(Fin[idx], fin(I+idx-ghosts))
+ exprs += (e,)
+ for idx in np.ndindex(*gr):
+ e = Assignment(fout(gr*I+idx), Fout_values[idx])
exprs += (e,)
- ekg.elementwise_kernel('interpolate_grid', *exprs,
- compute_resolution=dFin.compute_resolution, debug=True)
- raise NotImplementedError
+
+ interpolate_grid_kernel, _ = ekg.elementwise_kernel('interpolate_grid',
+ *exprs, compute_resolution=dFin.compute_resolution, debug=False)
+ exchange_ghosts = self.dFout.exchange_ghosts(build_launcher=True)
+
+ kl = OpenClKernelListLauncher(name='lowpass_filter')
+ kl += interpolate_grid_kernel
+ kl += exchange_ghosts
+
+ self.execute_kernels = functools.partial(kl, queue=self.cl_env.default_queue)
@op_apply
def apply(self, **kwds):
super(OpenClPolynomialInterpolationFilter, self).apply(**kwds)
- raise NotImplementedError
+ evt = self.execute_kernels()
class OpenClSubgridRestrictionFilter(SubgridRestrictionFilterBase, OpenClSymbolic):
diff --git a/hysop/numerics/interpolation/polynomial.py b/hysop/numerics/interpolation/polynomial.py
index a57a2973a..8aa98499a 100644
--- a/hysop/numerics/interpolation/polynomial.py
+++ b/hysop/numerics/interpolation/polynomial.py
@@ -96,7 +96,7 @@ class PolynomialInterpolationMethod(SpaceDiscretizationMethod):
Operator helper to handle polynomial interpolation method.
"""
__default_method = {
- MultiScaleInterpolation: PolynomialInterpolation.CUBIC_FDC6,
+ MultiScaleInterpolation: PolynomialInterpolation.NONIC_FDC2,
}
__available_methods = {
@@ -467,6 +467,8 @@ class PolynomialSubgridInterpolator(object):
----------
dim: int
Number of dimensions to interpolate (same as interpolator.dim).
+ ghosts: int
+ Number of required ghosts.
n: int
Corresponds to 2*(ghosts+1), the number of required nodes to generate the
polynomial coefficients (same as interpolator.n).
@@ -519,6 +521,7 @@ class PolynomialSubgridInterpolator(object):
P = interpolator.P
n = interpolator.n
N = interpolator.N
+ ghosts = interpolator.ghosts
s = tuple(gr+1 for gr in grid_ratio)
S = np.prod(s, dtype=np.int64)
@@ -531,6 +534,7 @@ class PolynomialSubgridInterpolator(object):
self.S = S
self.n = n
self.N = N
+ self.ghosts = interpolator.ghosts
self.dim = dim
self.interpolator = interpolator
self.key = key
diff --git a/hysop/symbolic/array.py b/hysop/symbolic/array.py
index 9e99a8fe6..4c28c0765 100644
--- a/hysop/symbolic/array.py
+++ b/hysop/symbolic/array.py
@@ -192,7 +192,7 @@ class SymbolicBuffer(SymbolicMemoryObject):
class SymbolicNdBuffer(SymbolicBuffer):
"""Same as a SymbolicBuffer, but with ndindex access and ghost support."""
- def __new__(cls, name, memory_object=None, dim=None, strides=None, ghosts=None, **kwds):
+ def __new__(cls, name, memory_object=None, dim=None, strides=None, dtype=None, ghosts=None, **kwds):
obj = super(SymbolicNdBuffer, cls).__new__(cls, memory_object=memory_object,
name=name, **kwds)
msg='Dimension could not be deduced from memory_object, '
@@ -216,21 +216,25 @@ class SymbolicNdBuffer(SymbolicBuffer):
pretty_name='g'+subscript(i),
dtype=npw.int32) for i in xrange(dim))
obj._allow_update_symbolic_constants = True
- obj.update_symbolic_constants(memory_object=memory_object, strides=strides, ghosts=ghosts, force=False)
+ obj.update_symbolic_constants(memory_object=memory_object, strides=strides, dtype=dtype, ghosts=ghosts, force=False)
return obj
- def bind_memory_object(self, memory_object, strides=None, ghosts=None, force=False, **kwds):
+ def bind_memory_object(self, memory_object, strides=None, dtype=None, ghosts=None, force=False, **kwds):
super(SymbolicNdBuffer, self).bind_memory_object(memory_object=memory_object, force=force, **kwds)
- self.update_symbolic_constants(memory_object=memory_object, strides=strides, ghosts=ghosts, force=force)
+ self.update_symbolic_constants(memory_object=memory_object, strides=strides, dtype=dtype, ghosts=ghosts, force=force)
- def update_symbolic_constants(self, memory_object, strides, ghosts, force):
+ def update_symbolic_constants(self, memory_object, strides, dtype, ghosts, force):
if hasattr(self, '_allow_update_symbolic_constants') and self._allow_update_symbolic_constants:
- strides = first_not_None(strides, getattr(memory_object, 'strides', None))
+ strides = first_not_None(strides, getattr(memory_object, 'strides', None))
+ dtype = first_not_None(dtype, getattr(memory_object, 'dtype', None))
assert (strides is not None), 'Could not determine strides from memory_object.'
+ assert (dtype is not None), 'Could not determine dtype from memory_object.'
+ itemsize = dtype.itemsize
strides = to_tuple(strides)
check_instance(strides, tuple, values=(int,long), size=self._dim)
for ss,si in zip(self._symbolic_strides, strides):
- ss.bind_value(si, force=force)
+ assert si%itemsize == 0
+ ss.bind_value(si//itemsize, force=force)
ghosts = first_not_None(ghosts, getattr(memory_object, 'ghosts', None))
assert (ghosts is not None), 'Could not determine ghosts from memory_object.'
diff --git a/hysop/tools/sympy_utils.py b/hysop/tools/sympy_utils.py
index 6e4634f3c..d656de1e4 100644
--- a/hysop/tools/sympy_utils.py
+++ b/hysop/tools/sympy_utils.py
@@ -154,7 +154,7 @@ class AppliedUndef(sm.function.AppliedUndef):
#return '{}({})'.format(self._pretty_name,
#','.join(printer._print(a) for a in self.args))
-def subscript(i, with_sign=False):
+def subscript(i, with_sign=False, disable_unicode=False):
"""
Generate an unicode subscript of value i, signs can be enforced.
"""
@@ -164,16 +164,19 @@ def subscript(i, with_sign=False):
s0 = snumber[0]
if s0 in decimals:
snumber = '+'+snumber
- out = u''
- for s in snumber:
- if s in decimals:
- out += decimal_subscripts[int(s)]
- elif s=='+':
- out += signs[0]
- elif s=='-':
- out += signs[1]
- else:
- out += s
+ if disable_unicode:
+ out = snumber
+ else:
+ out =u''
+ for s in snumber:
+ if s in decimals:
+ out += decimal_subscripts[int(s)]
+ elif s=='+':
+ out += signs[0]
+ elif s=='-':
+ out += signs[1]
+ else:
+ out += s
return out
def exponent(i, with_sign=False):
@@ -198,16 +201,21 @@ def exponent(i, with_sign=False):
out += s
return out
-def subscripts(ids,sep,with_sign=False,with_parenthesis=False,prefix=''):
+def subscripts(ids,sep,with_sign=False,with_parenthesis=False,prefix='',disable_unicode=False):
"""
Generate a unicode tuple subscript separated by sep,
with or without parenthesis, prefix, and signs.
"""
ids = to_tuple(ids)
if with_parenthesis:
- return u'{}{}{}{}'.format(prefix,parenthesis[0],sep.join([subscript(i,with_sign) for i in ids]),parenthesis[1])
+ lparen = '(' if disable_unicode else parenthesis[0]
+ rparen = ')' if disable_unicode else parenthesis[1]
+ base = '{}{}{}{}' if disable_unicode else u'{}{}{}{}'
+ return base.format(prefix,lparen,sep.join([subscript(i,with_sign,disable_unicode)
+ for i in ids]),rparen)
else:
- return u'{}{}'.format(prefix,sep.join([subscript(i,with_sign) for i in ids]))
+ base = '{}{}' if disable_unicode else u'{}{}'
+ return base.format(prefix,sep.join([subscript(i,with_sign,disable_unicode) for i in ids]))
def exponents(ids,sep,with_sign=False,with_parenthesis=False,prefix=''):
"""
--
GitLab