diff --git a/hysop/backend/device/opencl/opencl_array.py b/hysop/backend/device/opencl/opencl_array.py
index bb50347f6afd089181302a745d8d9d25fbb6da5c..caae51deaece578600a1fb8d128fd1755024d499 100644
--- a/hysop/backend/device/opencl/opencl_array.py
+++ b/hysop/backend/device/opencl/opencl_array.py
@@ -39,9 +39,19 @@ class OpenClArray(Array):
backend.check_queue(handle.queue)
super(OpenClArray,self).__init__(handle=handle, backend=backend, **kargs)
- def s(self, name, **kwds):
+ def as_symbolic_array(self, name, **kwds):
+ """
+ Return a symbolic array variable that contain a reference to this array.
+ """
from hysop.symbolic.array import OpenClSymbolicArray
- return OpenClSymbolicArray(array=self, name=name, **kwds)
+ return OpenClSymbolicArray(memory_object=self, name=name, **kwds)
+
+ def as_symbolic_buffer(self, name, **kwds):
+ """
+ Return a symbolic buffer variable that contain a reference to this array.
+ """
+ from hysop.symbolic.array import OpenClSymbolicBuffer
+ return OpenClSymbolicBuffer(memory_object=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 577847c1876e133c8c947b7b964ad7aec4d87075..394e748943f81ef4f759a0b0ff7e7a2e9b99650b 100644
--- a/hysop/backend/device/opencl/opencl_symbolic.py
+++ b/hysop/backend/device/opencl/opencl_symbolic.py
@@ -16,13 +16,16 @@ from hysop.tools.numpywrappers import npw
from hysop.fields.continuous_field import Field
from hysop.parameters.tensor_parameter import TensorParameter
from hysop.operator.base.custom_symbolic_operator import ValidExpressions, \
- SymbolicExpressionParser, CustomSymbolicOperatorBase
+ SymbolicExpressionParser, CustomSymbolicOperatorBase
from hysop.backend.device.opencl.opencl_operator import OpenClOperator
-from hysop.constants import ComputeGranularity, SpaceDiscretization, TranspositionState, DirectionLabels, SymbolicExpressionKind
+from hysop.constants import ComputeGranularity, SpaceDiscretization, TranspositionState, \
+ DirectionLabels, SymbolicExpressionKind
from hysop.numerics.interpolation.interpolation import Interpolation
from hysop.numerics.odesolvers.runge_kutta import TimeIntegrator
from hysop.topology.cartesian_descriptor import CartesianTopologyDescriptors
-from hysop.backend.device.opencl.autotunable_kernels.custom_symbolic import OpenClAutotunableCustomSymbolicKernel
+from hysop.backend.device.opencl.autotunable_kernels.custom_symbolic import \
+ OpenClAutotunableCustomSymbolicKernel
+from hysop.tools.sympy_utils import subscript, subscripts
class OpenClSymbolic(OpenClOperator):
"""
@@ -30,7 +33,8 @@ class OpenClSymbolic(OpenClOperator):
that require custom symbolic kernels.
"""
__default_method = CustomSymbolicOperatorBase._CustomSymbolicOperatorBase__default_method
- __available_methods = CustomSymbolicOperatorBase._CustomSymbolicOperatorBase__available_methods
+ __available_methods = \
+ CustomSymbolicOperatorBase._CustomSymbolicOperatorBase__available_methods
@classmethod
def default_method(cls):
@@ -44,6 +48,7 @@ class OpenClSymbolic(OpenClOperator):
am.update(cls.__available_methods)
return am
+
@debug
def __init__(self, **kwds):
"""
@@ -71,7 +76,56 @@ class OpenClSymbolic(OpenClOperator):
check_instance(name, str)
check_instance(exprs, tuple, values=ValidExpressions, minsize=1)
self.expressions[name] = exprs
+
+ def symbolic_buffers(self, *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))
+ 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
+ varnames = names
+ latexnames = names
+ for (name, var_name, latex_name) in zip(snames, var_names, latex_names):
+ check_instance(name, (str,unicode))
+ buf = OpenClSymbolicBuffer(name=name, var_name=var_name, latex_name=latex_name,
+ memory_object=None, **kwds)
+ buffers += (buf,)
+ return buffers
+
+ def symbolic_arrays(self, *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
+ varnames = names
+ latexnames = names
+ for (name, var_name, latex_name) in zip(snames, var_names, latex_names):
+ check_instance(name, (str,unicode))
+ arr = OpenClSymbolicArray(name=name, var_name=var_name, latex_name=latex_name,
+ memory_object=None, **kwds)
+ arrays += (arr,)
+ return arrays
+ def symbolic_tmp_scalars(self, *names, **kwds):
+ from hysop.symbolic.tmp import TmpScalar
+ scalars = ()
+ for name in names:
+ check_instance(name, str)
+ tmp = TmpScalar(name=name, **kwds)
+ scalars += (tmp,)
+ return scalars
+
@debug
def handle_method(self, method):
"""
@@ -133,7 +187,7 @@ class OpenClSymbolic(OpenClOperator):
_cmp('output_fields', output_fields, expr_info.output_fields, exprs)
_cmp('input_params', input_params, expr_info.input_params, exprs)
_cmp('output_params', output_params, expr_info.output_params, exprs)
- assert (0 <= self.cr <= expr_info.max_granularity), cr
+ assert (0 <= self.cr <= expr_info.max_granularity), self.cr
expr_info.compute_granularity = self.cr
expr_info.time_integrator = self.time_integrator
expr_info.interpolation = self.interpolation
@@ -149,7 +203,7 @@ class OpenClSymbolic(OpenClOperator):
for expr_info in self.expr_infos.values():
expr_info.extract_field_requirements()
- dim = expr_info.domain.dim
+ dim = expr_info.dim
field_reqs = expr_info.field_requirements
direction = expr_info.direction
has_direction = expr_info.has_direction
@@ -162,11 +216,15 @@ class OpenClSymbolic(OpenClOperator):
min_ghosts_per_components = {}
- for (fields, is_input, iter_requirements) in zip((self.input_fields, self.output_fields), (True, False),
- (requirements.iter_input_requirements, requirements.iter_output_requirements)):
+ for (fields, is_input, iter_requirements) in \
+ zip((self.input_fields, self.output_fields), (True, False),
+ (requirements.iter_input_requirements,
+ requirements.iter_output_requirements)):
if not fields:
continue
for (field, td, req) in iter_requirements():
+ if (field not in expr_info.input_fields):
+ continue
min_ghosts = npw.int_zeros(shape=(field.nb_components, field.dim))
if has_direction:
req.axes = axes
@@ -190,7 +248,8 @@ class OpenClSymbolic(OpenClOperator):
req.max_ghosts = npw.minimum(_max_ghosts, req.max_ghosts)
min_ghosts[index] = _min_ghosts.copy()
else:
- req.min_ghosts[dim-1-direction] = max(G, req.min_ghosts[dim-1-direction])
+ req.min_ghosts[dim-1-direction] = \
+ max(G, req.min_ghosts[dim-1-direction])
min_ghosts[index][dim-1-direction] = G
assert req.min_ghosts[dim-1-direction] >= G
assert req.max_ghosts[dim-1-direction] >= G
@@ -200,7 +259,8 @@ class OpenClSymbolic(OpenClOperator):
expr_info.min_ghosts = {k:npw.max(v, axis=0) for (k,v) in \
min_ghosts_per_components.iteritems()}
- expr_info.min_ghosts_per_components = {field:gpc[:,-1-direction] for (field,gpc) in min_ghosts_per_components.iteritems()}
+ expr_info.min_ghosts_per_components = {field:gpc[:,-1-direction]
+ for (field,gpc) in min_ghosts_per_components.iteritems()}
return requirements
@debug
@@ -231,7 +291,8 @@ class OpenClSymbolic(OpenClOperator):
build_opts=build_opts, autotuner_config=autotuner_config)
for (name, expr_info) in self.expr_infos.iteritems():
- kernel, args_dict, update_input_parameters = kernel_autotuner.autotune(expr_info=expr_info)
+ kernel, args_dict, update_input_parameters = \
+ kernel_autotuner.autotune(expr_info=expr_info)
kl = kernel.build_launcher(**args_dict)
self.symbolic_kernels[name] = (kl, update_input_parameters)
diff --git a/hysop/backend/device/opencl/operator/poisson_rotational.py b/hysop/backend/device/opencl/operator/poisson_rotational.py
index f71551ff6a22538d94fce1252a78cf785a431b74..d5ee401c5ef81f0a4e6553a499d70de75a253d59 100644
--- a/hysop/backend/device/opencl/operator/poisson_rotational.py
+++ b/hysop/backend/device/opencl/operator/poisson_rotational.py
@@ -20,20 +20,20 @@ class OpenClPoissonRotational(PoissonRotationalOperatorBase, OpenClSymbolic):
def initialize(self, **kwds):
# request the projection kernel if required
- if (projection != FieldProjection.NONE):
+ if (self.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')
+ I = local_indexes_symbols[:3]
+ K = self.symbolic_buffers('K', count=3)
+ W = self.symbolic_arrays('W', count=3)
+ divW, = self.symbolic_tmp_scalars('divW')
- expr = Assignment(divW, sum(K[i][I[i]]*W[i]))
+ expr = Assignment(divW, sum(K[i][I[i]]*W[i] for i in range(3)))
exprs = (expr,)
for i in xrange(3):
- expr = Assignment(W[i], W[i] - K[i]*diwW)
+ expr = Assignment(W[i], W[i] - K[i][I[i]]*divW)
exprs += (expr,)
-
- self.require_symbolic_kernel('solenoidal_reprojection_W', exprs)
+
+ self.require_symbolic_kernel('solenoidal_reprojection_W', *exprs)
self.reprojection_buffers = (K,W)
@@ -91,7 +91,7 @@ class OpenClPoissonRotational(PoissonRotationalOperatorBase, OpenClSymbolic):
K2d[start:end] = k2
start = end
if (self.projection != FieldProjection.NONE):
- self.reprojection_buffers[0][i].bind_buffer(K1d[start:end])
+ self.reprojection_buffers[0][i].bind_memory_object(K1d[start:end])
self.K1d = K1d
self.K2d = K2d
@@ -120,16 +120,20 @@ class OpenClPoissonRotational(PoissonRotationalOperatorBase, OpenClSymbolic):
@debug
def setup(self, work):
+ fft_buffer, = work.get_buffer(self, 'R2C_C2R')
+ if (self.projection != FieldProjection.NONE):
+ N = npw.prod(sizes, dtype=npw.int64)
+ for i in xrange(3):
+ self.reprojection_buffers[1][i].bind_memory_object(fft_buffer[i*N:(i+1)*N])
super(OpenClPoissonRotational, self).setup(work)
- self._build_fft_plans(work)
+ self._build_fft_plans(fft_buffer)
self._build_ghost_exchanger()
- def _build_fft_plans(self, work):
+ def _build_fft_plans(self, fft_buffer):
axes = self.axes
context = self.backend.cl_env.context
queue = self.backend.cl_env.default_queue
- fft_buffer = work.get_buffer(self, 'R2C_C2R')[0]
forward_plans, backward_plans = [],[]
for (i,Wi) in enumerate(self.W_buffers):
diff --git a/hysop/backend/host/host_array.py b/hysop/backend/host/host_array.py
index c8f03287d93fc332a893a2902b46171ec0a013b9..d333f137e2faa29a00492a075669be080b19e17f 100644
--- a/hysop/backend/host/host_array.py
+++ b/hysop/backend/host/host_array.py
@@ -50,9 +50,19 @@ class HostArray(Array):
# array interface
self.__array_interface__ = handle.__array_interface__
- def s(self, name, **kwds):
+ def as_symbolic_array(self, name, **kwds):
+ """
+ Return a symbolic array variable that contain a reference to this array.
+ """
from hysop.symbolic.array import HostSymbolicArray
- return HostSymbolicArray(array=self, name=name, **kwds)
+ return HostSymbolicArray(memory_object=self, name=name, **kwds)
+
+ def as_symbolic_buffer(self, name, **kwds):
+ """
+ Return a symbolic buffer variable that contain a reference to this array.
+ """
+ from hysop.symbolic.array import HostSymbolicBuffer
+ return HostSymbolicBuffer(memory_object=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 7795758896b4eea7332ea93a40b5208d1931c021..39912e4486072b744bb5423af042a147625a2cb7 100644
--- a/hysop/core/arrays/array.py
+++ b/hysop/core/arrays/array.py
@@ -136,12 +136,19 @@ class Array(object):
return self.backend._call(f, *args, **kargs)
@abstractmethod
- def s(self, name, **kwds):
+ def as_symbolic_array(self, name, **kwds):
"""
Return a symbolic array variable that contain a reference to this array.
"""
pass
+ @abstractmethod
+ def as_symbolic_buffer(self, name, **kwds):
+ """
+ Return a symbolic buffer variable that contain a reference to this array.
+ """
+ pass
+
@abstractmethod
@required_property
def get_ndim(self):
diff --git a/hysop/operator/base/custom_symbolic_operator.py b/hysop/operator/base/custom_symbolic_operator.py
index 57b0c6e59929b8fa44be1fa5907509210afcad79..2afebee94d44885863220765bf5b99c89a1d12a8 100644
--- a/hysop/operator/base/custom_symbolic_operator.py
+++ b/hysop/operator/base/custom_symbolic_operator.py
@@ -17,11 +17,14 @@ from hysop.topology.cartesian_descriptor import CartesianTopologyDescriptors
from hysop.symbolic import Dummy, space_symbols, dspace_symbols, time_symbol
from hysop.symbolic.relational import Assignment
from hysop.symbolic.misc import TimeIntegrate
+from hysop.symbolic.tmp import TmpScalar
+from hysop.symbolic.array import SymbolicArray, SymbolicBuffer
from hysop.symbolic.field import AppliedSymbolicField, SymbolicDiscreteField
from hysop.symbolic.parameter import SymbolicScalarParameter
from hysop.symbolic.misc import ApplyStencil
-from hysop.constants import ComputeGranularity, SpaceDiscretization, TranspositionState, DirectionLabels, SymbolicExpressionKind
+from hysop.constants import ComputeGranularity, SpaceDiscretization, TranspositionState, \
+ DirectionLabels, SymbolicExpressionKind
from hysop.numerics.odesolvers.runge_kutta import TimeIntegrator, ExplicitRungeKutta, Euler, RK2, RK3, RK4
from hysop.numerics.interpolation.interpolation import Interpolation
from hysop.numerics.stencil.stencil_generator import StencilGenerator, CenteredStencilGenerator, MPQ
@@ -34,14 +37,30 @@ class ExprDiscretizationInfo(object):
self.write_counter = {}
self.parameters = {}
- def read(self, field, index, count=1):
+ def read(self, obj, index, count=1):
check_instance(count, int)
- check_instance(field, DiscreteFieldView)
- self.read_counter.setdefault(field.dfield, npw.int_zeros(shape=(field.nb_components,)))[index] += count
- def write(self, field, index, count=1):
+ if isinstance(obj, DiscreteFieldView):
+ self.read_counter.setdefault(obj.dfield,
+ npw.int_zeros(shape=(obj.nb_components,)))[index] += count
+ elif isinstance(obj, (SymbolicArray, TmpScalar)):
+ self.read_counter.setdefault(obj, 0)
+ self.read_counter[obj] += 1
+ else:
+ msg='Unsupported type {}.'.format(type(obj))
+ raise TypeError(msg)
+
+ def write(self, obj, index=None, count=1):
check_instance(count, int)
- check_instance(field, DiscreteFieldView)
- self.write_counter.setdefault(field.dfield, npw.int_zeros(shape=(field.nb_components,)))[index] += count
+ if isinstance(obj, DiscreteFieldView):
+ self.write_counter.setdefault(obj.dfield,
+ npw.int_zeros(shape=(obj.nb_components,)))[index] += count
+ elif isinstance(obj, (SymbolicArray, TmpScalar)):
+ self.write_counter.setdefault(obj, 0)
+ self.write_counter[obj] += 1
+ else:
+ msg='Unsupported type {}.'.format(type(obj))
+ raise TypeError(msg)
+
def copy(self):
obj = ExprDiscretizationInfo()
for (field, counts) in self.read_counter.iteritems():
@@ -52,15 +71,25 @@ class ExprDiscretizationInfo(object):
counter += counts
obj.push_parameters(**self.parameters)
return obj
+
def update(self, other):
check_instance(other, ExprDiscretizationInfo)
- for (field, counts) in other.read_counter.iteritems():
- counter = self.read_counter.setdefault(field, npw.int_zeros(shape=(field.nb_components,)))
- counter += counts
- for (field, counts) in other.write_counter.iteritems():
- counter = self.write_counter.setdefault(field, npw.int_zeros(shape=(field.nb_components,)))
- counter += counts
+ for (obj, counts) in other.read_counter.iteritems():
+ if isinstance(obj, AppliedSymbolicField):
+ counter = self.read_counter.setdefault(field, npw.int_zeros(shape=(field.nb_components,)))
+ counter += counts
+ else:
+ self.read_counter.setdefault(obj, 0)
+ self.read_counter[obj] += counts
+ for (obj, counts) in other.write_counter.iteritems():
+ if isinstance(obj, AppliedSymbolicField):
+ counter = self.write_counter.setdefault(obj, npw.int_zeros(shape=(field.nb_components,)))
+ counter += counts
+ else:
+ self.write_counter.setdefault(obj, 0)
+ self.write_counter[obj] += counts
self.push_parameters(**other.parameters)
+
def push_parameters(self, *param, **kwd_params):
self.parameters.update(**kwd_params)
for p in param:
@@ -97,10 +126,13 @@ class SymbolicExpressionInfo(object):
# continuous part
self.domain = None
+ self.input_arrays = {}
+ self.output_arrays = {}
+ self.input_buffers = {}
self.input_fields = {}
self.output_fields = {}
- self.input_params = {}
- self.output_params = {}
+ self.input_params = {}
+ self.output_params = {}
self.direction = None
self.has_direction = None
@@ -125,18 +157,18 @@ class SymbolicExpressionInfo(object):
self.input_dfields = None
self.output_dfields = None
self.inout_dfields = None
- self.max_granularity = None
self.discretization_info = None
self.stencils = None
self.tmp_vars = None
+ self._dim = None
def _is_discretized(self):
"""Return true if the SymbolicExpressionInfo was discretized."""
return (self.dexprs is not None)
def _get_dim(self):
"""Shortcut to domain dimension."""
- assert (self.domain is not None)
- return self.domain.dim
+ assert (self._dim is not None)
+ return self._dim
def _get_fields(self):
"""Return input and output fields."""
fields = {k: v for (k,v) in self.input_fields.iteritems()}
@@ -147,6 +179,10 @@ class SymbolicExpressionInfo(object):
fields = {k: v for (k,v) in self.input_params.iteritems()}
fields.update(self.output_params)
return fields
+
+ @property
+ def max_granularity(self):
+ return (self.dim - 1)
dim = property(_get_dim)
fields = property(_get_fields)
@@ -159,17 +195,41 @@ class SymbolicExpressionInfo(object):
def check_field(self, field):
"""
Check if given continuous field is compatible with previously
- parsed fields.
+ parsed fields and arrays.
"""
check_instance(field, Field)
if (self.domain is None):
self.domain = field.domain
- self.max_granularity = self.dim - 1
- else:
- if (field.domain is not self.domain):
- msg='Domain mismatch for field {}:\n{}\nReference domain was:\n{}.'
- msg=msg.format(field.name, field.domain, self.domain)
+ if (self._dim is None):
+ self._dim = field.domain.dim
+ elif (self._dim != field.domain.dim):
+ msg='Dimension mismatch between field domain dimension and array dimension, '
+ msg+='got {} and {}.'.format(self._dim, field.domain.dim)
+ raise ValueError(msg)
+ elif (field.domain is not self.domain):
+ msg='Domain mismatch for field {}:\n{}\nReference domain was:\n{}.'
+ msg=msg.format(field.name, field.domain, self.domain)
+ raise ValueError(msg)
+
+ def check_array(self, array):
+ """
+ Check if given symbolic array is compatible with previously
+ parsed fields and arrays.
+ """
+ check_instance(array, SymbolicArray)
+ dim = array.dim
+ if (self.domain is not None):
+ if (self.domain.dim != dim):
+ msg='Dimension mismatch between field domain dimension and array dimension, '
+ msg+='got {} and {}.'.format(self.domain.dim, dim)
+ raise ValueError(msg)
+ elif (self._dim is not None):
+ if (self._dim != dim):
+ msg='Dimension mismatch between arrays, got {} and {}.'
+ msg=msg.format(self._dim, dim)
raise ValueError(msg)
+ else:
+ self._dim = dim
def extract_field_requirements(self):
field_requirements = SymbolicExpressionParser.extract_field_requirements(self)
@@ -180,12 +240,12 @@ class SymbolicExpressionInfo(object):
def discretize_expressions(self, input_dfields, output_dfields):
check_instance(input_dfields, dict, keys=Field, values=DiscreteFieldView)
check_instance(output_dfields, dict, keys=Field, values=DiscreteFieldView)
- assert set(input_dfields.keys()) == set(self.input_fields.keys())
- assert set(output_dfields.keys()) == set(self.output_fields.keys())
- self.input_dfields = input_dfields
- self.output_dfields = output_dfields
+ assert len(set(self.input_fields.keys()) - set(input_dfields.keys())) == 0
+ assert len(set(self.output_fields.keys()) - set(output_dfields.keys())) == 0
+ self.input_dfields = {k:v for (k,v) in input_dfields.iteritems() if (k in self.input_fields)}
+ self.output_dfields = {k:v for (k,v) in output_dfields.iteritems() if (k in self.output_fields)}
self.inout_dfields = {k:v for (k,v) in self.output_dfields.iteritems()
- if ((k in input_dfields) and (input_dfields[k].dfield is v.dfield))}
+ if ((k in self.input_dfields) and (self.input_dfields[k].dfield is v.dfield))}
self.stencils = {}
SymbolicExpressionParser.discretize_expressions(self)
@@ -215,7 +275,8 @@ class SymbolicExpressionInfo(object):
xd = space_symbols[direction]
if (xd != var):
msg='Expression already contained a derivative with respect to {} (direction {}).'
- msg+='\nFound a new derivative direction which is not compatible with the current one.'
+ msg+='\nFound a new derivative direction which is not compatible with the current '
+ msg+='one.'
msg+='\nCannot differentiate with respect to {}.'
msg=msg.format(xd, direction, var)
raise ValueError(msg)
@@ -238,8 +299,8 @@ class SymbolicExpressionInfo(object):
'''.format(
self.kind,
'\n'+'\n'.join(' {}/ {}'.format(i,e) for i,e in enumerate(self.exprs)),
- ', '.join('{}'.format(f.name) for f in self.input_fields.keys()) if self.input_fields else 'none',
- ', '.join('{}'.format(f.name) for f in self.output_fields.keys()) if self.output_fields else 'none',
+ ', '.join('{}'.format(f.name) for f in self.input_fields.keys()) if self.input_fields else 'none',
+ ', '.join('{}'.format(f.name) for f in self.output_fields.keys()) if self.output_fields else 'none',
', '.join('{}'.format(p) for p in self.input_params.keys()) if self.input_params else 'none',
', '.join('{}'.format(p) for p in self.output_params.keys()) if self.output_params else 'none',
'\n'+'\n'.join(' {}: {}'.format(f.name, (d.short_description() if isinstance(d, TopologyView)
@@ -249,8 +310,10 @@ class SymbolicExpressionInfo(object):
min_ghosts_per_components:{}
min_ghosts:{}
'''.format(
- '\n'+'\n'.join(' {}/ [{}]'.format(f.name,', '.join(str(x) for x in gpc)) for f,gpc in self.min_ghosts_per_components.iteritems()),
- '\n'+'\n'.join(' {}/ [{}]'.format(f.name,', '.join(str(x) for x in g)) for f,g in self.min_ghosts.iteritems()))
+ '\n'+'\n'.join(' {}/ [{}]'.format(f.name,', '.join(str(x) for x in gpc))
+ for f,gpc in self.min_ghosts_per_components.iteritems()),
+ '\n'+'\n'.join(' {}/ [{}]'.format(f.name,', '.join(str(x) for x in g))
+ for f,g in self.min_ghosts.iteritems()))
if self.is_discretized:
msg += '''
@@ -259,8 +322,10 @@ class SymbolicExpressionInfo(object):
write_counter: {}
'''.format(
'\n'+'\n'.join(' {}/ {}'.format(i,e) for i,e in enumerate(self.dexprs)),
- ', '.join('{}{}'.format(f.name, self.discretization_info.read_counter[f.dfield]) for f in self.input_dfields.values()) if self.input_dfields else 'none',
- ', '.join('{}{}'.format(f.name, self.discretization_info.write_counter[f.dfield]) for f in self.output_dfields.values()) if self.output_dfields else 'none')
+ ', '.join('{}{}'.format(f.name, self.discretization_info.read_counter[f.dfield])
+ for f in self.input_dfields.values()) if self.input_dfields else 'none',
+ ', '.join('{}{}'.format(f.name, self.discretization_info.write_counter[f.dfield])
+ for f in self.output_dfields.values()) if self.output_dfields else 'none')
return msg
@@ -275,11 +340,15 @@ class SymbolicExpressionParser(object):
def check_expressions(cls, exprs):
kind = None
fields = {}
+ arrays = {}
for expr in exprs:
check_instance(expr, ValidExpressions)
lhs = expr.args[0]
- field = lhs
- if isinstance(lhs, AppliedSymbolicField):
+ field = None
+ array = None
+ if isinstance(lhs, TmpScalar):
+ continue
+ elif isinstance(lhs, (AppliedSymbolicField, SymbolicArray)):
if (kind is None):
kind = SymbolicExpressionKind.AFFECT
elif (kind is not SymbolicExpressionKind.AFFECT):
@@ -287,17 +356,24 @@ class SymbolicExpressionParser(object):
msg+=' an expression that is of kind {}.\n expr: {}'
msg=msg.format(kind, SymbolicExpressionKind.AFFECT, expr)
raise ValueError(msg)
+ if isinstance(lhs, AppliedSymbolicField):
+ field = lhs
+ else:
+ array = lhs
elif isinstance(lhs, sm.Derivative):
- field = lhs.args[0]
- _vars = lhs.args[1:]
- _t = time_symbol
+ lhs, _vars = lhs.args[0], lhs.args[1:]
+ _t = time_symbol
unique_vars = set(_vars)
- if not isinstance(field, AppliedSymbolicField):
+ if isinstance(lhs, AppliedSymbolicField):
+ field = lhs
+ elif isinstance(lhs, SymbolicArray):
+ array = lhs
+ else:
msg='Assignment LHS cannot be a derivative of a {}.'
msg=msg.format(type(field))
raise TypeError(msg)
if (_t not in unique_vars) or len(unique_vars)!=1:
- msg='Assignment LHS only be a derivative of time {}, got {}.'
+ msg='Assignment LHS can only be a derivative of time {}, got {}.'
msg=msg.format(_t, ', '.join(str(x) for x in unique_vars))
raise TypeError(msg)
if (kind is None):
@@ -310,16 +386,31 @@ class SymbolicExpressionParser(object):
else:
msg='Assignment LHS cannot be of type {}.'.format(type(lhs))
raise TypeError(msg)
+
+ if (field is not None):
+ assert isinstance(field, AppliedSymbolicField)
+ index = field.index
+ field = field.field
+ key = (field, index)
+ if (key in fields):
+ msg='Field {} was already written by expression\n'
+ msg+='{}\ncannot write it again in expression\n'
+ msg+='{}\nFATAL ERROR: Invalid expressions.'
+ msg=msg.format(field.name, fields[key], expr)
+ raise ValueError(msg)
+ fields[key] = expr
+
+ if (array is not None):
+ assert isinstance(array, SymbolicArray)
+ key = array
+ if (key in arrays):
+ msg='Zrray {} was already written by expression\n'
+ msg+='{}\ncannot write it again in expression\n'
+ msg+='{}\nFATAL ERROR: Invalid expressions.'
+ msg=msg.format(array.name, arrays[key], expr)
+ raise ValueError(msg)
+ arrays[key] = expr
- assert isinstance(field, AppliedSymbolicField)
- index = field.index
- field = field.field
- key = (field, index)
- if (key in fields):
- msg='Field {} was already written by expression\n {}\ncannot write it again in expression\n {}\nFATAL ERROR: Invalid expressions.'
- msg=msg.format(field.name, fields[key], expr)
- raise ValueError(msg)
- fields[key] = expr
return kind
@classmethod
@@ -327,6 +418,16 @@ class SymbolicExpressionParser(object):
info = SymbolicExpressionInfo(name, exprs, **kwds)
for expr in exprs:
cls.parse_one(variables, info, expr)
+ if (info._dim is None):
+ msg='\n\nFATAL ERROR: Neither SymbolicFields nor SymbolicArrays were present in parsed '
+ msg+='symbolic expressions.'
+ msg+='\nAt least one is needed to deduce the shape of the compute kernel.'
+ msg+='\n'
+ msg+='\nExpressions were:'
+ for (i,e) in enumerate(exprs):
+ msg+='\n {:2>}/ {}'.format(i,e)
+ msg+='\n'
+ raise RuntimeError(msg)
if (info.direction is None):
info.direction = 0
return info
@@ -343,12 +444,11 @@ class SymbolicExpressionParser(object):
@classmethod
def parse_assignment(cls, variables, info, lhs, rhs):
- if isinstance(lhs, AppliedSymbolicField):
+ if isinstance(lhs, (AppliedSymbolicField, SymbolicArray, TmpScalar)):
cls.write(variables, info, lhs)
cls.parse_subexpr(variables, info, rhs)
elif isinstance(lhs, sm.Derivative):
f = lhs.args[0]
- assert isinstance(f, AppliedSymbolicField), f
cls.read(variables, info, f)
cls.write(variables, info, f)
cls.parse_subexpr(variables, info, rhs)
@@ -379,7 +479,16 @@ class SymbolicExpressionParser(object):
@classmethod
def write(cls, variables, info, var):
- if isinstance(var, AppliedSymbolicField):
+ if isinstance(var, TmpScalar):
+ return
+ elif isinstance(var, SymbolicArray):
+ array = var
+ info.check_array(array)
+ if (array not in info.output_arrays):
+ info.output_arrays[array.name] = array
+ else:
+ assert info.output_arrays[array.name] is array
+ elif isinstance(var, AppliedSymbolicField):
field = var.field
info.check_field(field)
if (field not in variables):
@@ -441,9 +550,12 @@ class SymbolicExpressionParser(object):
lhs = expr.args[0]
if isinstance(lhs, sm.Derivative):
lhs = lhs.args[0]
+ if not isinstance(lhs, AppliedSymbolicField):
+ continue
field, index = lhs.field, lhs.index
updated_fields[i] = (field, index, '{}_{}'.format(field.name, index))
- min_ghosts_per_expr[i] = {'{}_{}'.format(k[0].name,k[1]):v.min_ghosts[-direction-1] for (k,v) in field_reqs.iteritems()}
+ min_ghosts_per_expr[i] = {'{}_{}'.format(k[0].name,k[1]):v.min_ghosts[-direction-1]
+ for (k,v) in field_reqs.iteritems()}
for (field, reqs) in field_reqs.iteritems():
if field in field_requirements:
field_requirements[field].update_requirements(reqs)
@@ -474,29 +586,32 @@ class SymbolicExpressionParser(object):
nsteps = info.time_integrator.stages
else:
nsteps = 1
+
+ info.nsteps = nsteps
+ info.nlhsfields = nlhsfields
- G_f = expr_ghost_map[:,:nlhsfields]
- min_ghosts_per_step = npw.int_zeros(shape=(nsteps, nlhsfields))
- min_ghosts_per_step[0,:] = npw.max(G_f, axis=0)
- for s in xrange(1, nsteps):
- min_ghosts_per_step[s] = npw.max(G_f + min_ghosts_per_step[s-1][:,None], axis=0)
- min_ghosts_lhs = min_ghosts_per_step[nsteps-1]
-
- if nsteps>1:
- min_ghosts_rhs = npw.max(expr_ghost_map[:,nlhsfields:] + min_ghosts_per_step[nsteps-2][:,None], axis=0)
- else:
- min_ghosts_rhs = npw.max(expr_ghost_map[:,nlhsfields:], axis=0)
- min_ghosts = min_ghosts_lhs.tolist() + min_ghosts_rhs.tolist()
+ if (nlhsfields > 0):
+ G_f = expr_ghost_map[:,:nlhsfields]
+ min_ghosts_per_step = npw.int_zeros(shape=(nsteps, nlhsfields))
+ min_ghosts_per_step[0,:] = npw.max(G_f, axis=0)
+ for s in xrange(1, nsteps):
+ min_ghosts_per_step[s] = npw.max(G_f + min_ghosts_per_step[s-1][:,None], axis=0)
+ min_ghosts_lhs = min_ghosts_per_step[nsteps-1]
+
+ if nsteps>1:
+ min_ghosts_rhs = npw.max(expr_ghost_map[:,nlhsfields:] + \
+ min_ghosts_per_step[nsteps-2][:,None], axis=0)
+ else:
+ min_ghosts_rhs = npw.max(expr_ghost_map[:,nlhsfields:], axis=0)
+ min_ghosts = min_ghosts_lhs.tolist() + min_ghosts_rhs.tolist()
+ info.min_ghosts_per_integration_step = min_ghosts_per_step
+ info.min_ghosts_per_field_name = dict(zip(all_fields, min_ghosts))
+ info.expr_ghost_map = expr_ghost_map
lhs_fields = { v:k for (k,v) in lhs_fields.iteritems() }
lhs_fields = tuple(lhs_fields[i] for i in xrange(nlhsfields))
all_fields = { v:k for (k,v) in all_fields.iteritems() }
all_fields = tuple(all_fields[i] for i in xrange(nfields))
- info.nsteps = nsteps
- info.nlhsfields = nlhsfields
- info.expr_ghost_map = expr_ghost_map
- info.min_ghosts_per_integration_step = min_ghosts_per_step
- info.min_ghosts_per_field_name = dict(zip(all_fields, min_ghosts))
info.lhs_field_names = lhs_fields
info.rhs_field_names = rhs_fields
info.all_field_names = all_fields
@@ -529,7 +644,8 @@ class SymbolicExpressionParser(object):
if (direction is not None):
xd = space_symbols[direction]
if (unique_dvars - set([xd])):
- msg='Expression already contained a derivative with respect to {} (direction {}, {}-axis).'
+ msg='Expression already contained a derivative with respect to {} (direction {}, '
+ msg+='{}-axis).'
msg+='\nFound a new derivative direction which is not compatible with the current one.'
msg+='\nCannot differentiate with respect to {}.'
msg=msg.format(xd, direction, DirectionLabels[direction],
@@ -612,9 +728,10 @@ class SymbolicExpressionParser(object):
def discretize_assignment(cls, info, expr):
lhs, rhs = expr.args
rhs, di = cls.discretize_subexpr(info, rhs)
- if isinstance(lhs, AppliedSymbolicField):
+ if isinstance(lhs, (AppliedSymbolicField, SymbolicArray, TmpScalar,)):
func = expr.func
elif isinstance(lhs, sm.Derivative):
+ assert isinstance(lhs.args[0], AppliedSymbolicField)
assert lhs.args[1] == time_symbol
assert len(lhs.args)==2
lhs = lhs.args[0]
@@ -626,13 +743,16 @@ class SymbolicExpressionParser(object):
msg='Invalid symbolic assignment lhs type {}.'
msg=msg.format(type(lhs))
raise NotImplementedError(msg)
+
+ if isinstance(lhs, AppliedSymbolicField):
+ field, index, indexed_field = lhs.field, lhs.index, lhs.indexed_field
+ dfield = info.output_dfields[field]
+ lhs = dfield.s[index]
+ check_instance(lhs, SymbolicDiscreteField)
+ cls.write_discrete(info, lhs, dfield, di)
+ else:
+ di.write(lhs)
- field, index, indexed_field = lhs.field, lhs.index, lhs.indexed_field
- dfield = info.output_dfields[field]
- lhs = dfield.s[index]
- check_instance(lhs, SymbolicDiscreteField)
-
- cls.write_discrete(info, lhs, dfield, di)
new_expr = func(lhs, rhs)
return new_expr, di
@@ -801,7 +921,8 @@ class CustomSymbolicOperatorBase(DirectionalOperatorBase):
splitting_dim = first_not_None(splitting_dim, expr_info.domain.dim)
if (expr_info.direction != splitting_direction):
- msg='Direction mismatch, expression has derivative in direction {} but direction {} has been specified.'
+ msg='Direction mismatch, expression has derivative in direction {} but direction {} '
+ msg+='has been specified.'
msg=msg.format(expr_info.direction, direction)
raise ValueError(msg)
@@ -846,7 +967,8 @@ class CustomSymbolicOperatorBase(DirectionalOperatorBase):
min_ghosts_per_components = {}
- for (fields, is_input, iter_requirements) in zip((self.input_fields, self.output_fields), (True, False),
+ for (fields, is_input, iter_requirements) in zip((self.input_fields, self.output_fields),
+ (True, False),
(requirements.iter_input_requirements, requirements.iter_output_requirements)):
if not fields:
continue
@@ -883,7 +1005,8 @@ class CustomSymbolicOperatorBase(DirectionalOperatorBase):
min_ghosts_per_components[field] = min_ghosts
expr_info.min_ghosts = {k:npw.max(v, axis=0) for (k,v) in \
min_ghosts_per_components.iteritems()}
- expr_info.min_ghosts_per_components = {field:gpc[:,-1-direction] for (field,gpc) in min_ghosts_per_components.iteritems()}
+ expr_info.min_ghosts_per_components = {field:gpc[:,-1-direction]
+ for (field,gpc) in min_ghosts_per_components.iteritems()}
return requirements
diff --git a/hysop/operator/base/poisson_rotational.py b/hysop/operator/base/poisson_rotational.py
index f6b44f1d1fa07059546026f684723f4869cab848..9ef41ddc90bf1200f7e4b271de0c6310374ed540 100644
--- a/hysop/operator/base/poisson_rotational.py
+++ b/hysop/operator/base/poisson_rotational.py
@@ -86,6 +86,7 @@ class PoissonRotationalOperatorBase(FftOperatorBase):
self.W = vorticity
self.U = velocity
self.projection = projection
+ self.dim = velocity.dim
@debug
def discretize(self):
@@ -113,7 +114,6 @@ class PoissonRotationalOperatorBase(FftOperatorBase):
self.W_buffers = W_buffers
self.U_buffers = U_buffers
- self.dim = dim
self.dtype = dtype
self.ctype = ctype
self.axes = axes
diff --git a/hysop/symbolic/array.py b/hysop/symbolic/array.py
index 5425522878cd79b2fd5e8f96a686ffaf3596a96a..7ce473f75db6ca9b01671a57049b2cd7acf7b75e 100644
--- a/hysop/symbolic/array.py
+++ b/hysop/symbolic/array.py
@@ -7,43 +7,90 @@ 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)
+class SymbolicMemoryObject(SymbolicScalar):
+ def __new__(cls, memory_object, name, **kwds):
+ obj = super(SymbolicMemoryObject, cls).__new__(cls, name=name, **kwds)
+ obj._memory_object = None
+ if (memory_object is not None):
+ obj.bind_memory_object(memory_object)
return obj
+
+ @property
+ def is_bound(self):
+ return (self._memory_object is not None)
@property
- def array(self):
- if (self._array is None):
- msg='Symbolic array {} has not been setup yet.'.format(self.name)
+ def memory_object(self):
+ if (self._memory_object is None):
+ msg='Symbolic memory_object {} has not been setup yet.'.format(self.name)
raise RuntimeError(msg)
- return self._array
+ return self._memory_object
@abstractmethod
- def bind_array(self, array):
- if (self._array is not None):
- msg='An array was already bind to SymbolicArray {}.'.format(self.name)
+ def bind_memory_object(self, memory_object):
+ if (self._memory_object is not None):
+ msg='An memory_object was already bind to SymbolicArray {}.'.format(self.name)
raise RuntimeError(msg)
- if isinstance(array, (OpenClArray, HostArray)):
- array = array.handle
- self._array = array
+ if isinstance(memory_object, (OpenClArray, HostArray)):
+ memory_object = memory_object.handle
+ self._memory_object = memory_object
+
+class IndexedBuffer(sm.Indexed):
+ """
+ Tag for indexed SymbolicBuffers.
+ """
+ pass
+
+class SymbolicArray(SymbolicMemoryObject):
+ """
+ An array is considered to be indexed by local indices in
+ autogenerated vectorized code (in symbolic code generation
+ framework).
+ """
+ def __getitem__(self, key):
+ msg='Symbolic array {} cannot be indexed.'.format(self.name)
+ raise RuntimeError(msg)
+
+ @property
+ def array(self):
+ return self._memory_object
+
+ @property
+ def dim(self):
+ return 3
+class SymbolicBuffer(SymbolicMemoryObject):
+ """
+ A buffer will not be indexed by local indices by default.
+ The user has to index a SymbolicBuffer so that it can be
+ used for code generation.
+ """
def __getitem__(self, key):
assert isinstance(key, (int, npw.integer, sm.Expr))
- return sm.Indexed(self, key)
+ return IndexedBuffer(self, key)
+
+ @property
+ def buffer(self):
+ return self._memory_object
+
+class SymbolicHostMemoryObject(object):
+ def bind_memory_object(self, memory_object):
+ check_instance(memory_object, (HostArray, npw.ndarray))
+ super(SymbolicHostMemoryObject, self).bind_memory_object(memory_object)
-class HostSymbolicArray(SymbolicArray):
- def bind_array(self, array):
- check_instance(array, (HostArray, npw.ndarray))
- super(HostSymbolicArray, self).bind_array(array)
+class SymbolicDeviceMemoryObject(object):
+ def bind_memory_object(self, memory_object):
+ check_instance(memory_object, (OpenClArray, clArray.Array))
+ super(SymbolicDeviceMemoryObject, self).bind_memory_object(memory_object)
-class OpenClSymbolicArray(SymbolicArray):
- def bind_array(self, array):
- check_instance(array, (OpenClArray, clArray.Array))
- super(OpenClSymbolicArray, self).bind_array(array)
+class HostSymbolicArray(SymbolicHostMemoryObject, SymbolicArray):
+ pass
+class OpenClSymbolicArray(SymbolicDeviceMemoryObject, SymbolicArray):
+ pass
+class HostSymbolicBuffer(SymbolicHostMemoryObject, SymbolicBuffer):
+ pass
+class OpenClSymbolicBuffer(SymbolicDeviceMemoryObject, SymbolicBuffer):
+ pass
if __name__ == '__main__':
@@ -51,11 +98,27 @@ if __name__ == '__main__':
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]
+
+ A = HostSymbolicArray(a, name='a')
+ B = b.as_symbolic_array('b')
+
+ C = HostSymbolicBuffer(a, name='c')
+ D = b.as_symbolic_buffer('d')
+
+ print 7*A + 9*B + 10*C - 4*D
+
+ assert A.array is a
+ assert B.array is b.handle
+ assert C.buffer is a
+ assert D.buffer is b.handle
+
+ try:
+ A[5]
+ assert False
+ except RuntimeError as e:
+ pass
+
+ print C[5]
+ print D[C]
+
+ print
diff --git a/hysop/symbolic/tmp.py b/hysop/symbolic/tmp.py
new file mode 100644
index 0000000000000000000000000000000000000000..cf9e4aae54bb6fd0a62bce0faccea6e7f2702683
--- /dev/null
+++ b/hysop/symbolic/tmp.py
@@ -0,0 +1,5 @@
+
+from hysop.symbolic.base import SymbolicScalar
+class TmpScalar(SymbolicScalar):
+ """Tag for temporary scalar variables."""
+ pass
diff --git a/hysop/tools/types.py b/hysop/tools/types.py
index eb32a880eaafcf4049f9e44e98757313e4cb0edc..1ce7e9df84e9217fd6ac417b9cc72508198fb065 100644
--- a/hysop/tools/types.py
+++ b/hysop/tools/types.py
@@ -101,7 +101,8 @@ def check_instance(val, cls, allow_none=False, **kargs):
maxval = kargs.pop('maxval', None)
for v in val:
if not any((isinstance(v, val_cls) for val_cls in all_val_cls)):
- msg='Value contained in given {} is not an instance of {} but a value of type {}.'
+ msg='Value contained in given {} is not an instance of {} but a '
+ msg='value of type {}.'
msg=msg.format(type(val).__name__, all_val_cls, type(v))
print_offending_value(v, all_val_cls)
raise TypeError(msg)