diff --git a/hysop/backend/device/opencl/opencl_array.py b/hysop/backend/device/opencl/opencl_array.py
index 369c9f096f683d69b5150a2f2e66c35468ccec63..bb50347f6afd089181302a745d8d9d25fbb6da5c 100644
--- a/hysop/backend/device/opencl/opencl_array.py
+++ b/hysop/backend/device/opencl/opencl_array.py
@@ -38,6 +38,10 @@ class OpenClArray(Array):
backend.check_queue(handle.queue)
super(OpenClArray,self).__init__(handle=handle, backend=backend, **kargs)
+
+ def s(self, name, **kwds):
+ from hysop.symbolic.array import OpenClSymbolicArray
+ return OpenClSymbolicArray(array=self, name=name, **kwds)
def get_ndim(self):
return self._handle.ndim
diff --git a/hysop/backend/device/opencl/opencl_symbolic.py b/hysop/backend/device/opencl/opencl_symbolic.py
index 976ca37bf27a06e132f6f2d14351cd285fb5175e..577847c1876e133c8c947b7b964ad7aec4d87075 100644
--- a/hysop/backend/device/opencl/opencl_symbolic.py
+++ b/hysop/backend/device/opencl/opencl_symbolic.py
@@ -82,7 +82,7 @@ class OpenClSymbolic(OpenClOperator):
self.cr = method.pop(ComputeGranularity)
self.space_discretization = method.pop(SpaceDiscretization)
self.time_integrator = method.pop(TimeIntegrator)
- self.interpolation = method.pop(Interpolation)
+ self.interpolation = method.pop(Interpolation)
assert (2 <= self.space_discretization), self.space_discretization
assert (self.space_discretization % 2 == 0), self.space_discretization
diff --git a/hysop/backend/device/opencl/operator/poisson_rotational.py b/hysop/backend/device/opencl/operator/poisson_rotational.py
index ef9bdf36fe322a37b60aba5e94ea2f6fdc8b56c5..f71551ff6a22538d94fce1252a78cf785a431b74 100644
--- a/hysop/backend/device/opencl/operator/poisson_rotational.py
+++ b/hysop/backend/device/opencl/operator/poisson_rotational.py
@@ -9,12 +9,37 @@ from hysop.core.graph.graph import op_apply
from hysop.core.memory.memory_request import MemoryRequest
from hysop.backend.device.opencl.opencl_fft import OpenClFFT
from hysop.backend.device.codegen.base.variables import dtype_to_ctype
+from hysop.constants import FieldProjection
+from hysop.symbolic import local_indexes_symbols
+from hysop.symbolic.relational import Assignment
class OpenClPoissonRotational(PoissonRotationalOperatorBase, OpenClSymbolic):
'''
Solves the poisson-rotational equation using clFFT.
'''
+ def initialize(self, **kwds):
+ # request the projection kernel if required
+ if (projection != FieldProjection.NONE):
+ assert (self.dim == 3), self.dim
+ I = local_indexes_symbols[:3]
+ K = self.symbolic_buffers('Kx', 'Ky', 'Kz')
+ W = self.symbolic_arrays('Wx', 'Wy', 'Wz')
+ divW = self.tmp_scalar('divW')
+
+ expr = Assignment(divW, sum(K[i][I[i]]*W[i]))
+ exprs = (expr,)
+ for i in xrange(3):
+ expr = Assignment(W[i], W[i] - K[i]*diwW)
+ exprs += (expr,)
+
+ self.require_symbolic_kernel('solenoidal_reprojection_W', exprs)
+ self.reprojection_buffers = (K,W)
+
+
+
+ super(OpenClPoissonRotational, self).initialize(**kwds)
+
def generate_wave_numbers(self, dim, resolution, length, dtype, ctype, axes):
if (dim>3):
msg='clFFT only support 1D, 2D or 3D plans, got a {}D domain.'
@@ -57,7 +82,7 @@ class OpenClPoissonRotational(PoissonRotationalOperatorBase, OpenClSymbolic):
Ksizes, Koffsets = (), ()
start, end = 0,0
- for (k1,k2) in zip(K1,K2):
+ for i,(k1,k2) in enumerate(zip(K1,K2)):
assert k1.size == k2.size
Koffsets += (start,)
Ksizes += (k1.size,)
@@ -65,6 +90,8 @@ class OpenClPoissonRotational(PoissonRotationalOperatorBase, OpenClSymbolic):
K1d[start:end] = k1
K2d[start:end] = k2
start = end
+ if (self.projection != FieldProjection.NONE):
+ self.reprojection_buffers[0][i].bind_buffer(K1d[start:end])
self.K1d = K1d
self.K2d = K2d
diff --git a/hysop/backend/host/host_array.py b/hysop/backend/host/host_array.py
index 42be0548e1bfb810638163af2155548c4983bbe9..c8f03287d93fc332a893a2902b46171ec0a013b9 100644
--- a/hysop/backend/host/host_array.py
+++ b/hysop/backend/host/host_array.py
@@ -49,6 +49,10 @@ class HostArray(Array):
# array interface
self.__array_interface__ = handle.__array_interface__
+
+ def s(self, name, **kwds):
+ from hysop.symbolic.array import HostSymbolicArray
+ return HostSymbolicArray(array=self, name=name, **kwds)
def get_ndim(self):
return self.handle.ndim
diff --git a/hysop/core/arrays/array.py b/hysop/core/arrays/array.py
index 661cc6eb38f5ee9f8415424c4b6eca1e380c7674..7795758896b4eea7332ea93a40b5208d1931c021 100644
--- a/hysop/core/arrays/array.py
+++ b/hysop/core/arrays/array.py
@@ -134,7 +134,13 @@ class Array(object):
"""
f = getattr(self.handle, fname)
return self.backend._call(f, *args, **kargs)
-
+
+ @abstractmethod
+ def s(self, name, **kwds):
+ """
+ Return a symbolic array variable that contain a reference to this array.
+ """
+ pass
@abstractmethod
@required_property
diff --git a/hysop/operator/base/poisson_rotational.py b/hysop/operator/base/poisson_rotational.py
index b1825697a7675e83acc72a38757e60a249253afb..f6b44f1d1fa07059546026f684723f4869cab848 100644
--- a/hysop/operator/base/poisson_rotational.py
+++ b/hysop/operator/base/poisson_rotational.py
@@ -49,7 +49,8 @@ class PoissonRotationalOperatorBase(FftOperatorBase):
input_fields = { vorticity: wtopology }
output_fields = { velocity: vtopology }
- if (projection == FieldProjection.NONE) or (projection is None) or (velocity.dim==2):
+ if (projection == FieldProjection.NONE) or (projection is None) \
+ or (projection==0) or (velocity.dim==2):
projection = FieldProjection.NONE
output_fields[vorticity] = wtopology
else:
@@ -78,8 +79,6 @@ class PoissonRotationalOperatorBase(FftOperatorBase):
freq = projection
if (freq>=1):
self._do_project = lambda simu: ((simu.current_iteration % freq)==0)
- elif (freq==0):
- self._do_project = lambda simu: False
else:
msg='Got negative reprojection frequency {}.'.format(freq)
raise ValueError(msg)
diff --git a/hysop/symbolic/__init__.py b/hysop/symbolic/__init__.py
index 4dfab950398b174c75a55c1f67fb07b66437e6f6..bb113c31ad19f739bca9fc68faac526ff5ad1776 100644
--- a/hysop/symbolic/__init__.py
+++ b/hysop/symbolic/__init__.py
@@ -27,3 +27,9 @@ dspace_symbols = tuple(SpaceSymbol(name=u'dx'+subscript(i),
latex_name='dx_{{{}}}'.format(i))
for i in xrange(16))
"""Dummy symbols representing space infinitesimals."""
+
+local_indexes_symbols = tuple(SpaceSymbol(name=u'i'+subscript(i),
+ var_name='i{}'.format(i),
+ latex_name='i_{{{}}}'.format(i))
+ for i in xrange(16))
+"""Dummy symbols local array indexes."""
diff --git a/hysop/symbolic/array.py b/hysop/symbolic/array.py
new file mode 100644
index 0000000000000000000000000000000000000000..5425522878cd79b2fd5e8f96a686ffaf3596a96a
--- /dev/null
+++ b/hysop/symbolic/array.py
@@ -0,0 +1,61 @@
+
+from abc import ABCMeta, abstractmethod
+from hysop.symbolic.base import SymbolicScalar, sm
+from hysop.tools.types import check_instance, to_tuple, first_not_None
+from hysop.tools.numpywrappers import npw
+from hysop.backend.device.opencl import clArray
+from hysop.backend.device.opencl.opencl_array import OpenClArray
+from hysop.backend.host.host_array import HostArray
+
+class SymbolicArray(SymbolicScalar):
+ def __new__(cls, array, name, **kwds):
+ obj = super(SymbolicArray, cls).__new__(cls, name=name, **kwds)
+ obj._array = None
+ if (array is not None):
+ obj.bind_array(array)
+ return obj
+
+ @property
+ def array(self):
+ if (self._array is None):
+ msg='Symbolic array {} has not been setup yet.'.format(self.name)
+ raise RuntimeError(msg)
+ return self._array
+
+ @abstractmethod
+ def bind_array(self, array):
+ if (self._array is not None):
+ msg='An array was already bind to SymbolicArray {}.'.format(self.name)
+ raise RuntimeError(msg)
+ if isinstance(array, (OpenClArray, HostArray)):
+ array = array.handle
+ self._array = array
+
+ def __getitem__(self, key):
+ assert isinstance(key, (int, npw.integer, sm.Expr))
+ return sm.Indexed(self, key)
+
+class HostSymbolicArray(SymbolicArray):
+ def bind_array(self, array):
+ check_instance(array, (HostArray, npw.ndarray))
+ super(HostSymbolicArray, self).bind_array(array)
+
+class OpenClSymbolicArray(SymbolicArray):
+ def bind_array(self, array):
+ check_instance(array, (OpenClArray, clArray.Array))
+ super(OpenClSymbolicArray, self).bind_array(array)
+
+
+if __name__ == '__main__':
+ from hysop.deps import sm
+ from hysop.core.arrays.all import default_host_array_backend
+ a = npw.ones(shape=(10,10), dtype=npw.int8)
+ b = default_host_array_backend.zeros(shape=(10,), dtype=npw.uint16)
+ A = HostSymbolicArray(array=a, name='a')
+ B = b.s('b')
+ print 7*A + 9*B
+ assert A.array is a
+ assert B.array is b.handle
+
+ print A[5]
+ print A[B]