diff --git a/examples/particles_above_salt/particles_above_salt_symmetrized.py b/examples/particles_above_salt/particles_above_salt_symmetrized.py
index d26e0873dd427db6dca885594177e58325923864..621badee6c55a87720463f659ffffefcf3418397 100644
--- a/examples/particles_above_salt/particles_above_salt_symmetrized.py
+++ b/examples/particles_above_salt/particles_above_salt_symmetrized.py
@@ -30,7 +30,7 @@ def delta(Ys, l0):
def init_concentration(data, coords, l0):
X = coords[-1].copy()
- X = 600.0 - np.abs(X-1200.0)
+ X = np.abs(X-1200.0) - 600.0
Ys = coords[:-1]
data[0][...] = 0.5*(1.0 +
sp.special.erf((X-delta(Ys,l0))/l0))
@@ -132,12 +132,14 @@ def compute(args):
advected_fields = (vorti,S),
velocity_cfl = args.cfl,
variables = {velo: npts, vorti: npts, S: npts},
- dt=dt)
+ dt=dt, **extra_op_kwds)
+ # mirror sediment settling speed at Y=1200 (Y in [0..2400])
VP = [0]*dim
- VP[-1] = 'Select({},{},gidx.y<{})'.format(tg.dump(-Vp),
- tg.dump(+Vp),
- (npts[0]-1)//2)
+ VP[-1] = 'select({}, {}, field_infos->P.pos.X<{})'.format(
+ tg.dump(-Vp),
+ tg.dump(+Vp),
+ tg.dump(1200.0))
advec_C = DirectionalAdvection(implementation=impl,
name='advec_C',
@@ -146,7 +148,7 @@ def compute(args):
relative_velocity = VP,
velocity_cfl = args.cfl,
variables = {velo: npts, C: npts},
- dt=dt)
+ dt=dt, **extra_op_kwds)
#> Stretch and diffuse vorticity
if (dim==3):
@@ -348,8 +350,8 @@ if __name__=='__main__':
box_origin=(0.0,), box_length=(1.0,),
tstart=0.0, tend=500.0,
dt=1e-6, cfl=0.5, lcfl=0.125,
- dump_times=tuple(float(x) for x in range(500)),
- dump_freq=0)
+ dump_times=None,#tuple(float(x) for x in range(500)),
+ dump_freq=100)
parser.run(compute)
diff --git a/hysop/backend/device/codegen/base/variables.py b/hysop/backend/device/codegen/base/variables.py
index b96c6e458715eb6ae3fb0881cc2d17583bf35a59..3924583104befe1d7e844c9ec012d85006a77b7d 100644
--- a/hysop/backend/device/codegen/base/variables.py
+++ b/hysop/backend/device/codegen/base/variables.py
@@ -237,7 +237,7 @@ class CodegenVariable(object):
ptr_const = self.ptr_const if (ptr_const is None) else ptr_const
ptr_volatile = self.ptr_volatile if (ptr_volatile is None) else ptr_volatile
ptr_restrict = self.ptr_restrict if (ptr_restrict is None) else ptr_restrict
-
+
return cls(name=name, nl=nl,
value=value, svalue=svalue, init=init,
ctype=ctype, storage=storage,
@@ -914,6 +914,11 @@ class CodegenStruct(CodegenVariable):
struct_var=struct_var)
self.genvars(struct, var_overrides)
+
+
+ def newvar(self, *args, **kwds):
+ kwds.setdefault('cls', CodegenVariable)
+ return super(CodegenStruct, self).newvar(*args, **kwds)
def force_symbolic(self, force=True):
if force is None:
@@ -936,7 +941,7 @@ class CodegenStruct(CodegenVariable):
else:
raise AttributeError
- def genvars(self,struct,var_overrides):
+ def genvars(self, struct, var_overrides):
if var_overrides is None:
var_overrides = {}
diff --git a/hysop/backend/device/codegen/functions/advection_rhs.py b/hysop/backend/device/codegen/functions/advection_rhs.py
index 9a4e7fe5bc9981e3713098d024ca4103bb377da3..08ce04c49eb2556daae7b9839202977cabd3f8ca 100644
--- a/hysop/backend/device/codegen/functions/advection_rhs.py
+++ b/hysop/backend/device/codegen/functions/advection_rhs.py
@@ -21,7 +21,8 @@ class DirectionalAdvectionRhsFunction(OpenClFunctionCodeGenerator):
relative_velocity,
ptr_restrict=True,
itype='int',
- known_args=None):
+ known_args=None,
+ field_infos=None):
assert work_dim>=1 and work_dim<=3
check_instance(boundary,BoundaryCondition)
@@ -38,7 +39,7 @@ class DirectionalAdvectionRhsFunction(OpenClFunctionCodeGenerator):
vtype = typegen.vtype(ftype,nparticles)
(args,basename) = self.build_prototype(typegen,work_dim,itype,ftype,vtype,
- nparticles,ptr_restrict,storage,is_cached, is_periodic)
+ nparticles,ptr_restrict,storage,is_cached, is_periodic, field_infos)
reqs = self.build_requirements(typegen,work_dim,itype,ftype,vtype,
nparticles,ptr_restrict,storage,is_cached, is_periodic)
@@ -65,7 +66,7 @@ class DirectionalAdvectionRhsFunction(OpenClFunctionCodeGenerator):
self.gencode()
def build_prototype(self, typegen, work_dim, itype, ftype, vtype,
- nparticles, ptr_restrict, storage, is_cached, is_periodic):
+ nparticles, ptr_restrict, storage, is_cached, is_periodic, field_infos):
args = ArgDict()
@@ -74,8 +75,6 @@ class DirectionalAdvectionRhsFunction(OpenClFunctionCodeGenerator):
nl=True)
args['X'] = CodegenVectorClBuiltin('X', ftype, nparticles, typegen, add_impl_const=True,
nl=True)
- #args['gidx'] = CodegenVectorClBuiltin('gidx', ftype, nparticles, typegen, add_impl_const=True,
- #nl=True)
if is_periodic and (not is_cached):
args['line_width'] = CodegenVariable('line_width', itype, typegen,
@@ -87,6 +86,10 @@ class DirectionalAdvectionRhsFunction(OpenClFunctionCodeGenerator):
args['inv_dx'] = CodegenVariable('inv_dx', ftype, typegen, add_impl_const=True, nl=True)
# args['active'] = CodegenVariable('active','bool',typegen, add_impl_const=True)
+ if (field_infos is not None):
+ args['field_infos'] = field_infos.pointer(name='field_infos', ptr_level=1,
+ const=True, ptr_const=True, nl=True)
+
cached = 'cached_' if is_cached else ''
basename = 'advection_rhs_{}{}{}p'.format(cached, ftype[0], nparticles)
diff --git a/hysop/backend/device/codegen/kernels/directional_advection.py b/hysop/backend/device/codegen/kernels/directional_advection.py
index 75b17be9f5b7217d07c30643aa0310ad2282d200..5d66ff58af650d2e038bc4db54f469afc9296abc 100644
--- a/hysop/backend/device/codegen/kernels/directional_advection.py
+++ b/hysop/backend/device/codegen/kernels/directional_advection.py
@@ -24,8 +24,9 @@ from hysop.backend.device.opencl.opencl_array_backend import OpenClArrayBackend
from hysop.backend.device.codegen.base.utils import WriteOnceDict, ArgDict
from hysop.backend.device.codegen.base.statistics import WorkStatistics
-from hysop.backend.device.codegen.base.variables import CodegenStruct
-from hysop.backend.device.codegen.structs.mesh_info import MeshBaseStruct, MeshInfoStruct
+from hysop.backend.device.codegen.base.variables import CodegenStruct
+from hysop.backend.device.codegen.structs.mesh_info import MeshBaseStruct, MeshInfoStruct
+from hysop.backend.device.codegen.structs.indices import GlobalFieldInfos
from hysop.backend.device.codegen.functions.runge_kutta import RungeKuttaFunction
from hysop.backend.device.codegen.functions.advection_rhs import DirectionalAdvectionRhsFunction
@@ -194,13 +195,21 @@ class DirectionalAdvectionKernelGenerator(KernelCodeGenerator):
mesh_info_struct = MeshInfoStruct(typegen=typegen, vsize=vsize)
reqs['MeshInfoStruct'] = mesh_info_struct
+
+ field_names = ('V','P')
+ global_field_infos = GlobalFieldInfos(typegen=typegen, field_names=field_names,
+ workdim=work_dim, vsize=nparticles)
+ reqs['GlobalFieldInfos'] = global_field_infos
+ self.field_names = field_names
+ self.field_infos = global_field_infos.build_codegen_variable(name='field_infos')
+
advection_rhs = DirectionalAdvectionRhsFunction(typegen=typegen, work_dim=work_dim,
ftype=ftype, is_cached=is_cached,
boundary=vboundary[0], nparticles=nparticles,
relative_velocity=relative_velocity,
ptr_restrict=True,
- itype=itype)
+ itype=itype, field_infos=self.field_infos)
used_vars = RungeKuttaFunction._default_used_vars.copy()
used_vars['y']='X'
@@ -263,6 +272,7 @@ class DirectionalAdvectionKernelGenerator(KernelCodeGenerator):
storage = s.storage
nparticles = s.nparticles
min_ghosts = s.min_ghosts
+ field_infos = s.field_infos
symbolic_mode = s.symbolic_mode
use_short_circuit = s.use_short_circuit
@@ -271,7 +281,8 @@ class DirectionalAdvectionKernelGenerator(KernelCodeGenerator):
is_cached = s.is_cached
cache_size_known = s.cache_size_known
- vtype = tg.vtype(ftype,work_dim)
+ vtype = tg.vtype(ftype, work_dim)
+ pvtype = tg.vtype(ftype, nparticles)
global_id = s.vars['global_id']
local_id = s.vars['local_id']
@@ -307,13 +318,7 @@ class DirectionalAdvectionKernelGenerator(KernelCodeGenerator):
position_gid = CodegenVectorClBuiltin('P_gid', itype, work_dim, typegen=tg)
velocity_gid = CodegenVectorClBuiltin('V_gid', itype, work_dim, typegen=tg)
- global_ids = {}
- for fname in ['V', 'P']:
- global_ids[fname] = [CodegenVariable('{}_K{}'.format(fname,
- DirectionLabels[work_dim-1-i]),
- itype, const=True, typegen=tg) for i in xrange(work_dim-1) ]
- global_ids[fname] += [CodegenVectorClBuiltin('{}_KX'.format(fname), itype,
- nparticles, typegen=tg, const=True)]
+
runge_kutta = self.reqs['runge_kutta']
@@ -391,12 +396,16 @@ class DirectionalAdvectionKernelGenerator(KernelCodeGenerator):
line_work, advec_ghosts))
s.jumpline()
- for (k,v) in global_ids.iteritems():
+ for k in s.field_names:
+ v = field_infos[k]
+ v0 = '{}.idx.IX'.format(v)
+ v1 = '{}.pos.X'.format(v)
mi = s.vars['{}_mesh_info'.format(k)]
- v[-1].declare(codegen=s, init='{} + {} + {}'.format(
- mi['start'][0],
- line_offset,
- poffset))
+ s.append('{} = {} + {} + {};'.format(v0, mi['start'][0], line_offset, poffset))
+ s.append('{} = {} + convert_{}({})*{};'.format(v1, mi['global_mesh']['xmin'][0],
+ pvtype, v0, mi['dx'][0]))
+ s.jumpline()
+
position_gid.affect(i=0, codegen=s, init='{} + {}'.format(
line_offset, position_ghosts))
velocity_gid.affect(i=0, codegen=s, init='{} + {} - {}'.format(
@@ -415,10 +424,17 @@ class DirectionalAdvectionKernelGenerator(KernelCodeGenerator):
line_position.declare(al, init='{} + {}'.format(position,
line_position_offset), align=True)
else:
- for (k,v) in global_ids.iteritems():
+ for k in s.field_names:
+ v = field_infos[k]
+ d = DirectionLabels[i]
+ v0 = '{}.idx.I{}'.format(v, d)
+ v1 = '{}.pos.{}'.format(v, d)
mi = s.vars['{}_mesh_info'.format(k)]
- v[work_dim-1-i].declare(codegen=s,
- init='{} + {}'.format(mi['start'][i],'kji'[i]))
+ s.append('{} = {} + {};'.format(v0, mi['start'][i], 'kji'[i]))
+ s.append('{} = {} + {}*{};'.format(v1, mi['global_mesh']['xmin'][i],
+ v0, mi['dx'][i]))
+ s.jumpline()
+
s.append('{} = {} + {};'.format(position_gid[i], 'kji'[i],
position_ghosts))
s.append('{} = {} + {};'.format(velocity_gid[i], 'kji'[i],
@@ -429,6 +445,7 @@ class DirectionalAdvectionKernelGenerator(KernelCodeGenerator):
nested_loops = [_work_iterate_(i) for i in xrange(work_dim-1,-1,-1)]
with s._kernel_():
+ field_infos.declare(s)
s.jumpline()
with s._align_() as al:
global_id.declare(al,align=True,const=True)
@@ -491,7 +508,8 @@ class DirectionalAdvectionKernelGenerator(KernelCodeGenerator):
'X': X,
'dt': dt,
'inv_dx': inv_dx,
- 'line_offset': line_offset
+ 'line_offset': line_offset,
+ 'field_infos': '&{}'.format(field_infos),
}
if is_cached:
rk_args['line_velocity'] = Vc
diff --git a/hysop/backend/device/codegen/structs/indices.py b/hysop/backend/device/codegen/structs/indices.py
new file mode 100644
index 0000000000000000000000000000000000000000..dfe81ca3cd108b75189f6306b226ba23c41a02e9
--- /dev/null
+++ b/hysop/backend/device/codegen/structs/indices.py
@@ -0,0 +1,141 @@
+
+from hysop.deps import np
+from hysop.tools.types import check_instance, first_not_None
+from hysop.tools.misc import upper_pow2_or_3
+from hysop.backend.device.codegen.base.enum_codegen import EnumCodeGenerator
+from hysop.backend.device.codegen.base.struct_codegen import StructCodeGenerator
+from hysop.backend.device.opencl.opencl_types import OpenClTypeGen
+from hysop.constants import DirectionLabels
+
+class SpaceIndices(StructCodeGenerator):
+ def __init__(self, name, typegen, var_prefix, stype, workdim,
+ vsize=1, typedef=True, comments=None):
+ assert vsize in [1,2,3,4,8,16]
+ stype = stype.replace('fbtype', typegen.fbtype)
+ name += '{}D_{}{}'.format(workdim, stype.capitalize(), vsize)
+ if isinstance(typedef, bool):
+ if (typedef is True):
+ typedef = name+'_s'
+ else:
+ typedef = None
+ dtype = self.build_dtype(typegen, var_prefix, stype, vsize, workdim)
+ super(SpaceIndices, self).__init__(name=name,
+ dtype=dtype,
+ typegen=typegen,
+ typedef=typedef,
+ comments=comments)
+ self.workdim = workdim
+ self.vsize = vsize
+ self.stype = stype
+
+ @staticmethod
+ def build_dtype(typegen, var_prefix, stype, vsize, workdim):
+ tg = typegen
+ stype1 = tg.dtype_from_str(stype)
+ stypen = tg.dtype_from_str('{}{}'.format(stype, vsize))
+
+ dtype = []
+
+ fn = '{}{}'.format(var_prefix, DirectionLabels[0])
+ field = (fn, stypen)
+ dtype.append(field)
+
+ for i in xrange(workdim-2, -1, -1):
+ fn = '{}{}'.format(var_prefix, DirectionLabels[workdim-1-i])
+ field = (fn, stype1)
+ dtype.append(field)
+
+
+ return np.dtype(dtype)
+
+ def create(self, *args, **kwds):
+ raise NotImplementedError()
+
+
+class GlobalIndices(SpaceIndices):
+ def __init__(self, typegen, workdim, **kwds):
+ super(GlobalIndices, self).__init__(
+ typegen=typegen, workdim=workdim,
+ name='GlobalIndices', var_prefix='I', stype='int',
+ **kwds)
+
+class GlobalCoordinates(SpaceIndices):
+ def __init__(self, typegen, workdim, **kwds):
+ super(GlobalCoordinates, self).__init__(
+ typegen=typegen, workdim=workdim,
+ name='GlobalCoordinates', var_prefix='', stype='fbtype',
+ **kwds)
+
+class GlobalFieldInfo(StructCodeGenerator):
+ def __init__(self, typegen, workdim, vsize, typedef=True, **kwds):
+ gi = GlobalIndices(typegen=typegen, workdim=workdim, vsize=vsize)
+ gc = GlobalCoordinates(typegen=typegen, workdim=workdim, vsize=vsize)
+ dtype = self.build_dtype(gi, gc)
+ name='GlobalFieldInfo{}D_{}{}_{}{}'.format(workdim,
+ gi.stype.capitalize(), vsize,
+ gc.stype.capitalize(), vsize)
+ if isinstance(typedef, bool):
+ if (typedef is True):
+ typedef = name+'_s'
+ else:
+ typedef = None
+ super(GlobalFieldInfo, self).__init__(name=name,
+ dtype=dtype,
+ typedef=typedef,
+ typegen=typegen,
+ **kwds)
+ self.require(gi.name, gi)
+ self.require(gc.name, gc)
+ self.workdim = workdim
+ self.vsize = vsize
+
+ @staticmethod
+ def build_dtype(gi, gc):
+ dtype = [('idx', gi.dtype), ('pos', gc.dtype)]
+ return np.dtype(dtype)
+
+class GlobalFieldInfos(StructCodeGenerator):
+ def __init__(self, field_names, typegen, workdim, vsize, typedef=True, **kwds):
+ gfi = GlobalFieldInfo(typegen=typegen, workdim=workdim, vsize=vsize)
+ dtype = []
+ for fname in field_names:
+ field = (fname, gfi.dtype)
+ dtype.append(field)
+ dtype = np.dtype(dtype)
+ name = gfi.name.replace('GlobalFieldInfo', 'GlobalFieldInfos')
+ name += '_' + '_'.join(field_names)
+ if isinstance(typedef, bool):
+ if (typedef is True):
+ typedef = name+'_s'
+ else:
+ typedef = None
+ super(GlobalFieldInfos, self).__init__(name=name,
+ dtype=dtype,
+ typedef=typedef,
+ typegen=typegen,
+ **kwds)
+ self.require(gfi.name, gfi)
+ self.field_names = field_names
+ self.workdim = workdim
+ self.vsize = vsize
+
+
+if __name__ == '__main__':
+ from hysop.backend.device.codegen.base.opencl_codegen import OpenClCodeGenerator
+ from hysop.backend.device.codegen.base.test import _test_typegen
+
+ tg = _test_typegen('double', float_dump_mode='dec')
+
+ idx = GlobalIndices(typegen=tg, workdim=3, vsize=4)
+ pos = GlobalCoordinates(typegen=tg, workdim=3, vsize=4)
+ gfi = GlobalFieldInfo(typegen=tg, workdim=3, vsize=2)
+ gfis = GlobalFieldInfos(typegen=tg, workdim=2, vsize=8, field_names=('P','V'))
+
+ cg = OpenClCodeGenerator('test_generator',tg)
+ cg.require('idx', idx)
+ cg.require('pos', pos)
+ cg.require('gfi', gfi)
+ cg.require('gfis', gfis)
+
+ cg.edit()
+ cg.test_compile()