From 784600301a595a7114378c07d7b7c31a9fae2261 Mon Sep 17 00:00:00 2001
From: Jean-Baptiste Keck <Jean-Baptiste.Keck@imag.fr>
Date: Sun, 14 Apr 2019 16:06:38 +0200
Subject: [PATCH] broken advection
---
.../backend/device/codegen/base/variables.py | 10 +-
.../device/codegen/functions/polynomial.py | 2 +-
.../codegen/kernels/directional_advection.py | 478 ++++++++----------
.../autotunable_kernels/advection_dir.py | 52 +-
.../operator/directional/advection_dir.py | 7 +-
.../operator/directional/advection_dir.py | 5 +-
hysop/fields/cartesian_discrete_field.py | 16 +-
hysop/fields/continuous_field.py | 14 +-
hysop/numerics/interpolation/polynomial.py | 3 +-
hysop/operator/base/advection_dir.py | 55 +-
.../operator/tests/test_bilevel_advection.py | 30 +-
.../tests/test_directional_advection.py | 38 +-
hysop/tools/parameters.py | 2 +-
13 files changed, 348 insertions(+), 364 deletions(-)
diff --git a/hysop/backend/device/codegen/base/variables.py b/hysop/backend/device/codegen/base/variables.py
index e10575633..2e84e7fae 100644
--- a/hysop/backend/device/codegen/base/variables.py
+++ b/hysop/backend/device/codegen/base/variables.py
@@ -528,11 +528,11 @@ class CodegenArray(CodegenVariable):
_svalue = [None]*s0
for d in xrange(s0):
_name = '{}_{}'.format(name, d)
- _value = value[d]
- _svalue = svalue[d]
- val, sval = CodegenArray.initialize_rec(_name, typegen,
+ dvalue = value[d]
+ dsvalue = svalue[d]
+ val, sval = CodegenArray._initialize_rec(_name, typegen,
storage, ctype, const, volatile,
- _shape, _sshape, _value, _svalue,
+ _shape, _sshape, dvalue, dsvalue,
_ptr_level, _ptr_restrict, _ptr_const, _ptr_volatile,
symbolic_mode)
var = CodegenArray(name=_name, typegen=typegen,
@@ -1031,7 +1031,7 @@ class CodegenStruct(CodegenVariable):
if __name__ == '__main__':
- from hysop.backend.device.codegen.base.test import test_typegen
+ from hysop.backend.device.codegen.base.test import _test_typegen as test_typegen
tg = test_typegen('float')
print':: ARRAYS ::'
diff --git a/hysop/backend/device/codegen/functions/polynomial.py b/hysop/backend/device/codegen/functions/polynomial.py
index 1cb592148..78b1a689d 100644
--- a/hysop/backend/device/codegen/functions/polynomial.py
+++ b/hysop/backend/device/codegen/functions/polynomial.py
@@ -95,7 +95,7 @@ if __name__ == '__main__':
tg = _test_typegen('float')
pf = PolynomialFunction(tg,'float',4,
- [0,1,2,3,4,0,0,0,5,6,7,8,9,0], 'test_poly', True)
+ [0,1,2,3,4,0,0,0,5,6,7,8,9,0], 'test_poly', 'x', True)
pf.edit()
print pf.per_work_statistics()
diff --git a/hysop/backend/device/codegen/kernels/directional_advection.py b/hysop/backend/device/codegen/kernels/directional_advection.py
index 43bef696b..aba77d14c 100644
--- a/hysop/backend/device/codegen/kernels/directional_advection.py
+++ b/hysop/backend/device/codegen/kernels/directional_advection.py
@@ -1,4 +1,5 @@
import contextlib
+import itertools as it
from contextlib import contextmanager
from hysop.deps import math, operator, hashlib
@@ -32,6 +33,7 @@ from hysop.backend.device.codegen.functions.runge_kutta import RungeKuttaFunc
from hysop.backend.device.codegen.functions.advection_rhs import DirectionalAdvectionRhsFunction
from hysop.numerics.odesolvers.runge_kutta import ExplicitRungeKutta
+from hysop.numerics.interpolation.polynomial import PolynomialInterpolator, PolynomialSubgridInterpolator
from hysop.backend.device.opencl import cl, clCharacterize
from hysop.backend.device.opencl.opencl_env import OpenClEnvironment
@@ -43,41 +45,69 @@ from hysop.fields.discrete_field import DiscreteScalarFieldView
class DirectionalAdvectionKernelGenerator(KernelCodeGenerator):
@staticmethod
- def codegen_name(ftype, is_cached, rk_scheme, nparticles, min_ghosts,
- relative_velocity, bilevel, **kargs):
+ def codegen_name(ftype, is_cached, rk_scheme, nparticles, cache_ghosts,
+ relative_velocity, bilevel, interpolator, **kargs):
cache = ''
if is_cached:
- cache = 'cached'
- if bilevel is not None:
- cache += str(bilevel[-1])
- cache += '_'
- return 'directional_{}advection_{}_{}{}p_{}g__{}'.format(cache,rk_scheme.name(),
- ftype[0],nparticles,min_ghosts,abs(hash(relative_velocity)))
+ cache = 'cached_'
+ if (bilevel is not None):
+ bilevel=''
+ if interpolator:
+ bilevel += ['linear','cubic','quintic','septic','nonic','unknown'][min((interpolator.interpolator.deg[0]-1)//2, 5)]
+ bilevel += '_fdc{}'.format(interpolator.interpolator.fd[0])
+ else:
+ bilevel += 'legacy_linear'
+ bilevel+='_bilevel_'
+ else:
+ bilevel=''
+ return '{}directional_{}advection_{}_{}{}p_{}g__{}'.format(
+ bilevel, cache,
+ rk_scheme.name(),
+ ftype[0], nparticles, cache_ghosts,
+ abs(hash(relative_velocity)))
def __init__(self, typegen, work_dim, ftype,
is_cached, rk_scheme,
vboundary, nparticles,
relative_velocity,
+ velocity_cfl, grid_ratio,
offset_by_xmin = False,
- min_ghosts = 0,
use_short_circuit = None,
unroll_loops = None,
symbolic_mode = False,
debug_mode = False,
tuning_mode = False,
known_vars = None,
- is_bilevel = None):
-
+ polynomial_interpolator = None): # <= overrides legacy linear interpolator with custom interpolator for known grid ratios
assert work_dim>0 and work_dim<=3
assert nparticles in [1,2,4,8,16]
assert isinstance(relative_velocity, (str,float))
check_instance(vboundary[0],BoundaryCondition)
check_instance(vboundary[1],BoundaryCondition)
check_instance(rk_scheme, ExplicitRungeKutta)
-
- if (is_bilevel is not None) and (nparticles>1):
- msg='Bilevel support with multiple particles at a time has not been implemented yet.'
- raise NotImplementedError(msg)
+ check_instance(grid_ratio, np.ndarray, minval=1.0)
+ check_instance(polynomial_interpolator, PolynomialSubgridInterpolator, allow_none=True)
+
+ assert grid_ratio.dtype in (np.float32, np.float64)
+ grid_ratio = grid_ratio[::-1] # <= we pass from hysop representation ZYX to opencl XYZ
+ inverse_grid_ratio = 1.0/grid_ratio
+ integer_grid_ratio = np.ceil(grid_ratio).astype(np.int32)
+
+ is_bilevel = (grid_ratio>1.0).any()
+ velocity_cache_ghosts = self.advection_cache_ghosts(velocity_cfl, grid_ratio[0])
+
+ if is_bilevel and (is_cached is False):
+ msg='Bilevel advection only works if cache is enabled.'
+ raise RuntimeError(msg)
+
+ if (polynomial_interpolator is not None):
+ if is_bilevel:
+ print 'LOL'
+ assert (integer_grid_ratio == grid_ratio).all()
+ assert (integer_grid_ratio[::-1] == polynomial_interpolator.gr).all()
+ else:
+ msg='A polynomial_interpolator has been supplied but without specifying grid_ratio.'
+ raise RuntimeError(msg)
known_vars = known_vars or dict()
@@ -91,9 +121,9 @@ class DirectionalAdvectionKernelGenerator(KernelCodeGenerator):
is_periodic = (vboundary[0]==BoundaryCondition.PERIODIC \
and vboundary[1]==BoundaryCondition.PERIODIC)
- assert (is_periodic and not is_cached) or min_ghosts>0
+ assert (is_periodic and not is_cached) or velocity_cache_ghosts>0
- cache_size_known = ('local_size' in known_vars or is_bilevel is not None)
+ cache_size_known = ('local_size' in known_vars)
_global = OpenClCodeGenerator.default_keywords['global']
_local = OpenClCodeGenerator.default_keywords['local']
@@ -108,8 +138,8 @@ class DirectionalAdvectionKernelGenerator(KernelCodeGenerator):
itype = 'int'
name = DirectionalAdvectionKernelGenerator.codegen_name(
- ftype, is_cached, rk_scheme, nparticles, min_ghosts,
- relative_velocity, is_bilevel)
+ ftype, is_cached, rk_scheme, nparticles, velocity_cache_ghosts,
+ relative_velocity, is_bilevel, polynomial_interpolator)
kernel_reqs = self.build_requirements(
typegen, work_dim, itype, ftype, is_cached,
@@ -143,34 +173,42 @@ class DirectionalAdvectionKernelGenerator(KernelCodeGenerator):
self.is_periodic = is_periodic
self.is_cached = is_cached
self.is_bilevel = is_bilevel
- self.min_ghosts = min_ghosts
+ self.grid_ratio = grid_ratio
+ self.integer_grid_ratio = integer_grid_ratio
+ self.inverse_grid_ratio = inverse_grid_ratio
+ self.velocity_cache_ghosts = velocity_cache_ghosts
self.tuning_mode = tuning_mode
self.offset_by_xmin = offset_by_xmin
self.relative_velocity = relative_velocity
self.use_short_circuit = use_short_circuit
self.unroll_loops = unroll_loops
+ self.polynomial_interpolator = polynomial_interpolator
self.gencode()
+ self.edit()
@classmethod
- def advec_ghosts(cls, velocity_cfl):
- """Return the minimal numbers of ghosts required on the lasr axe of the velocity grid."""
- return int(1+math.floor(velocity_cfl))
+ def advection_ghosts(cls, velocity_cfl):
+ """Return the minimal numbers of ghosts required on the last axe of the velocity grid."""
+ assert velocity_cfl > 0.0, 'cfl <= 0.0'
+ return int(math.floor(velocity_cfl))
@classmethod
- def min_ghosts(cls, work_dim, velocity_cfl):
- """Return the minimal numbers of ghosts required on the velocity grid."""
- assert velocity_cfl > 0.0, 'cfl <= 0.0'
- ghosts = [0]*work_dim
- ghosts[-1] = cls.advec_ghosts(velocity_cfl)
- return np.asarray(ghosts, dtype=np.int32)
+ def advection_cache_ghosts(cls, velocity_cfl, grid_ratio):
+ assert isinstance(grid_ratio, float)
+ assert grid_ratio >= 1.0
+ linear_interpolation_ghosts = 1
+ return linear_interpolation_ghosts+int(math.floor(grid_ratio*velocity_cfl))
@classmethod
- def min_wg_size(cls, work_dim, velocity_cfl):
- """Return the minimum workgroup size."""
- ghosts = cls.min_ghosts(work_dim, velocity_cfl)
- return np.asarray(2*ghosts+1, dtype=np.int32).copy()
+ def min_velocity_ghosts(cls, work_dim, velocity_cfl, grid_ratio, periodicity):
+ """Return the minimal numbers of ghosts required on the velocity grid, includint interpolation ghosts."""
+ interpolation_ghosts = 1
+ ghosts = periodicity * (grid_ratio > 1.0)
+ ghosts += interpolation_ghosts
+ ghosts[-1] += cls.advection_ghosts(velocity_cfl)
+ return np.asarray(ghosts, dtype=np.int32)
def required_workgroup_cache_size(self, local_work_size):
"""Return a tuple of required (static,dynamic,total) cache bytes per workgroup."""
@@ -183,8 +221,8 @@ class DirectionalAdvectionKernelGenerator(KernelCodeGenerator):
sc,dc = 0,0
if is_cached:
- count = self.nparticles*local_work_size[0]+2*self.min_ghosts
- if 'local_size' in self.known_vars:
+ count = self.nparticles*local_work_size[0]+2*self.velocity_cache_ghosts
+ if ('local_size' in self.known_vars):
assert (self.known_vars['local_size'] == local_work_size[:work_dim]).all()
sc += count
else:
@@ -196,13 +234,6 @@ class DirectionalAdvectionKernelGenerator(KernelCodeGenerator):
return (sc,dc,tc)
- def required_workgroup_velocity_cache_size(self):
- ftype = self.ftype
- flt_bytes = self.typegen.FLT_BYTES[ftype]
- c = self.is_bilevel[-1]+2*self.min_ghosts
- c *= flt_bytes
- return c
-
def build_requirements(self,typegen,work_dim,itype,ftype,is_cached,rk_scheme,
vboundary,relative_velocity,nparticles,force_symbolic,storage,is_periodic,known_vars):
tg=typegen
@@ -216,20 +247,12 @@ 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, field_infos=self.field_infos)
+ itype=itype)
used_vars = RungeKuttaFunction._default_used_vars.copy()
used_vars['y']='X'
@@ -291,9 +314,10 @@ class DirectionalAdvectionKernelGenerator(KernelCodeGenerator):
vboundary = s.vboundary
storage = s.storage
nparticles = s.nparticles
- min_ghosts = s.min_ghosts
- field_infos = s.field_infos
+ velocity_cache_ghosts = s.velocity_cache_ghosts
tuning_mode = s.tuning_mode
+ is_bilevel = s.is_bilevel
+ polynomial_interpolator = s.polynomial_interpolator
symbolic_mode = s.symbolic_mode
use_short_circuit = s.use_short_circuit
@@ -322,44 +346,43 @@ class DirectionalAdvectionKernelGenerator(KernelCodeGenerator):
velocity, velocity_strides = s.velocity, s.velocity_strides
position, position_strides = s.position, s.position_strides
- has_bilevel = (self.is_bilevel is not None)
-
- if has_bilevel:
- assert is_cached and cache_size_known, \
- "In bilevel, velocity must be cached " + str(is_cached) + " " + str(cache_size_known)
-
- compute_grid_size = position_mesh_info['local_mesh']['compute_resolution'].view(
+ p_grid_size = position_mesh_info['local_mesh']['resolution'].view(
+ 'p_grid_size', slice(0, work_dim), const=True)
+ p_compute_grid_size = position_mesh_info['local_mesh']['compute_resolution'].view(
'p_compute_grid_size', slice(0, work_dim), const=True)
- if has_bilevel:
- v_grid_size = velocity_mesh_info['local_mesh']['resolution'].view(
+ p_grid_ghosts = position_mesh_info['ghosts'].view(
+ 'p_grid_ghosts', slice(0,work_dim), const=True)
+ p_dx = position_mesh_info['dx'].view('p_dx', slice(0,work_dim), const=True)
+ p_inv_dx = position_mesh_info['inv_dx'].view('p_inv_dx', slice(0,work_dim), const=True)
+
+ v_grid_size = velocity_mesh_info['local_mesh']['resolution'].view(
'v_grid_size', slice(0, work_dim), const=True)
-
- P_grid_ghosts = position_mesh_info['ghosts'].view(
- 'P_grid_ghosts', slice(0,work_dim), const=True)
- V_grid_ghosts = velocity_mesh_info['ghosts'].view(
- 'V_grid_ghosts', slice(0,work_dim), const=True)
-
- dx = position_mesh_info['dx'].view('p_dx', slice(0,work_dim), const=True)
- inv_dx = position_mesh_info['inv_dx'].view('p_inv_dx', slice(0,work_dim), const=True)
+ v_compute_grid_size = velocity_mesh_info['local_mesh']['compute_resolution'].view(
+ 'v_compute_grid_size', slice(0, work_dim), const=True)
+ v_grid_ghosts = velocity_mesh_info['ghosts'].view(
+ 'v_grid_ghosts', slice(0,work_dim), const=True)
v_dx = velocity_mesh_info['dx'].view('v_dx', slice(0,work_dim), const=True)
v_inv_dx = velocity_mesh_info['inv_dx'].view('v_inv_dx', slice(0,work_dim), const=True)
- xmin = CodegenVariable(name='xmin', ctype=ftype, typegen=tg, const=True,
+ p_xmin = CodegenVariable(name='p_xmin', ctype=ftype, typegen=tg, const=True,
init='{} + {}*{}'.format(position_mesh_info['local_mesh']['xmin'][0],
- P_grid_ghosts[0], dx[0]))
+ p_grid_ghosts[0], p_dx[0]))
+ v_xmin = CodegenVariable(name='v_xmin', ctype=ftype, typegen=tg, const=True,
+ init='{} + {}*{}'.format(velocity_mesh_info['local_mesh']['xmin'][0],
+ v_grid_ghosts[0], v_dx[0]))
position_gid = CodegenVectorClBuiltin('P_gid', itype, work_dim, typegen=tg)
velocity_gid = CodegenVectorClBuiltin('V_gid', itype, work_dim, typegen=tg)
-
- if has_bilevel:
- velocity_gid_pos = CodegenVectorClBuiltin('V_pos', ftype, work_dim, typegen=tg)
+ if is_bilevel:
+ velocity_p = CodegenVectorClBuiltin('V_pos', ftype, work_dim, typegen=tg)
velocity_h = CodegenVectorClBuiltin('Vh', ftype, work_dim, typegen=tg)
- velocity_ix = CodegenVariable(name='V_gid_x', ctype=itype, typegen=tg, const=True)
- line_offset_for_v = CodegenVariable(name='line_offset_for_v', ctype=itype, typegen=tg, const=False)
runge_kutta = self.reqs['runge_kutta']
-
- advec_ghosts = CodegenVariable('V_advec_ghosts', itype, typegen=tg, value=min_ghosts)
+
+ grid_ratio = CodegenVectorClBuiltin('grid_ratio', ftype, work_dim, typegen=tg, value=self.grid_ratio)
+ igrid_ratio = CodegenVectorClBuiltin('igrid_ratio', ftype, work_dim, typegen=tg, value=self.inverse_grid_ratio)
+ V_cache_ghosts = CodegenVariable('V_cache_ghosts', itype, typegen=tg, init=velocity_cache_ghosts)
+
line_offset = CodegenVariable(name='line_offset', ctype=itype, typegen=tg, const=True)
line_work = CodegenVariable(name='line_work', ctype=itype, typegen=tg, const=True)
@@ -379,16 +402,8 @@ class DirectionalAdvectionKernelGenerator(KernelCodeGenerator):
V_cache_width = CodegenVariable('V_cache_width', itype, typegen=tg,
const=True)
if cache_size_known:
- if has_bilevel:
- # if bilevel, velocity cache is the entire local line
- # WARNING : is_bilevel is HySoP velocity resolution (using transposition state)
- # whereas 'local_size' is a OpenCL size
- # if bilevel, velocity cache is the entire local line
- L = s.is_bilevel[-1]
- velocity_cache_full_length = True
- else:
- L = s.known_vars['local_size'][0]
- Vc_shape = (nparticles*L+2*min_ghosts,)
+ L = s.known_vars['local_size'][0]
+ Vc_shape = (nparticles*L+2*velocity_cache_ghosts,)
Vc = CodegenArray(name='Vc',dim=1,ctype=ftype,typegen=tg,
shape=Vc_shape, storage='__local')
del L
@@ -403,7 +418,16 @@ class DirectionalAdvectionKernelGenerator(KernelCodeGenerator):
kmax = CodegenVariable('kmax', itype, tg, const=True,
init='({Sx}+{npart}*{Lx}-1)/({npart}*{Lx})'.format(
- Sx=compute_grid_size[0], npart=npart(),Lx=local_size[0]))
+ Sx=p_compute_grid_size[0], npart=npart(),Lx=local_size[0]))
+
+ if polynomial_interpolator:
+ pi = polynomial_interpolator
+ weights = polynomial_interpolator.Wr.reshape(pi.gr+pi.n)
+ interpolation_weights = CodegenArray(name='interpolation_weights',ctype=ftype, dim=2*work_dim,
+ value=weights, typegen=tg, storage='__constant', const=True)
+ widx = CodegenVectorClBuiltin('widx', itype, work_dim, typegen=tg)
+ with self._codeblock_('global_scope_constants'):
+ interpolation_weights.declare(s)
@contextmanager
def _work_iterate_(i):
@@ -415,9 +439,9 @@ class DirectionalAdvectionKernelGenerator(KernelCodeGenerator):
else:
fval = global_id[i]
gsize = global_size[i]
- N = '{Sx}'.format(Sx=compute_grid_size[i])
- position_ghosts = P_grid_ghosts[i]
- velocity_ghosts = V_grid_ghosts[i]
+ N = '{Sx}'.format(Sx=p_compute_grid_size[i])
+ position_ghosts = p_grid_ghosts[i]
+ velocity_ghosts = v_grid_ghosts[i]
with s._for_('int {i}={fval}; {i}<{N}; {i}+={gsize}'.format(
i='kji'[i], fval=fval, gsize=gsize,N=N),
@@ -429,37 +453,25 @@ class DirectionalAdvectionKernelGenerator(KernelCodeGenerator):
is_active.declare(al, align=True,
init='({k}*{ls}+{lid} <= ({gs}+{npart}-1)/{npart})'.format(
npart=npart, lid=local_id[0], k='kji'[0],
- gs=compute_grid_size[0], ls=local_size[0]))
+ gs=p_compute_grid_size[0], ls=local_size[0]))
s.jumpline()
with s._align_() as al:
line_offset.declare(al, align=True, init='{}*{}*{}'.format(
'kji'[0], local_size[0], npart))
line_work.declare(al, align=True, init='({} ? {}-{} : {}*{})'.format(
- is_last, compute_grid_size[0], line_offset, local_size[0], npart))
- if is_cached and not has_bilevel:
+ is_last, p_compute_grid_size[0], line_offset, local_size[0], npart))
+ if is_cached:
al.jumpline()
V_cache_width.declare(al, align=True, init='{} + 2*{}'.format(
- line_work, advec_ghosts))
- s.jumpline()
-
- 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)]
- 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]))
+ line_work, V_cache_ghosts))
s.jumpline()
position_gid.affect(i=0, codegen=s, init='{} + {}'.format(
line_offset, position_ghosts))
- if not has_bilevel:
+ if not is_bilevel:
velocity_gid.affect(i=0, codegen=s, init='{} + {} - {}'.format(
- line_offset, velocity_ghosts, advec_ghosts))
+ line_offset, velocity_ghosts, V_cache_ghosts))
line_velocity_offset = ' + '.join('{}*{}'.format(velocity_gid[j], velocity_strides[j])
for j in xrange(work_dim-1, -1, -1))
@@ -468,74 +480,18 @@ class DirectionalAdvectionKernelGenerator(KernelCodeGenerator):
with s._align_() as al:
line_position.declare(al, init='{} + {}'.format(position, line_position_offset), align=True)
- if (s.is_bilevel is None):
+ if not is_bilevel:
line_velocity.declare(al, init='{} + {}'.format(velocity, line_velocity_offset), align=True)
else:
- 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)]
- 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))
- if has_bilevel:
- s.append('{} = ({} * {}) * {};'.format(
- velocity_gid_pos[i], 'kji'[i], dx[i], v_inv_dx[i]))
- s.append('{} = convert_int_rtn({});'.format(velocity_gid[i], velocity_gid_pos[i]))
- s.append('{} = {} - convert_{}({});'.format(
- velocity_h[i], velocity_gid_pos[i], part_ftype, velocity_gid[i]))
- if i==1:
- if work_dim==3:
- line_velocity_offset = '({}+{})*{}+({}+{})*{}+{}-{}'.format(
- velocity_gid[2],V_grid_ghosts[2], velocity_strides[2],
- velocity_gid[1],V_grid_ghosts[1], velocity_strides[1],
- V_grid_ghosts[0], advec_ghosts)
- elif work_dim==2:
- line_velocity_offset = '({}+{})*{}+{}-{}'.format(
- velocity_gid[1],V_grid_ghosts[1], velocity_strides[1],
- V_grid_ghosts[0], advec_ghosts)
- s.jumpline()
- with s._align_() as al:
- line_velocity.declare(al, init='{} + {}'.format(velocity, line_velocity_offset), align=True)
- line_offset_for_v.declare(al, init='0')
- s.jumpline()
- if velocity_cache_full_length:
- s.comment("Load velocity in cache with linear interpolation")
- with s._for_('int {i}={fval}; {i}<{N}; {i}+={gsize}'.format(
- i=velocity_ix, fval=local_id[0], gsize=local_size[0],N=Vc_shape[0]),
- unroll=False):
- is_first = True
- if work_dim==3:
- for ii,jj in [(ii,jj) for ii in range(2) for jj in range(2)]:
- s.append('{}[{}] {} ({}{})*({}{})*{};'.format(
- Vc, velocity_ix, " =" if is_first else "+=",
- "1.0-" if ii==0 else " ", velocity_h[2],
- "1.0-" if jj==0 else " ", velocity_h[1],
- s.vload(n=nparticles, ptr=line_velocity,
- offset='({})*{}+({})*{}+{}*{}'.format(
- ii, velocity_strides[2],
- jj, velocity_strides[1],
- velocity_ix, velocity_strides[0]))))
- is_first = False
- elif work_dim==2:
- for ii in range(2):
- s.append('{}[{}] {} ({}{})*{};'.format(
- Vc, velocity_ix, " =" if is_first else "+=",
- "1.0-" if ii==0 else " ", velocity_h[1],
- s.vload(n=nparticles,ptr=line_velocity,
- offset='({})*{}+{}*{}'.format(
- ii, velocity_strides[1],
- velocity_ix, velocity_strides[0]))))
- is_first = False
- else:
- raise RuntimeError("Bilevel advection 1D need further developments")
- s.barrier(_local=True)
- s.jumpline()
+ if is_bilevel:
+ s.append('{} = {} * {} + {};'.format(
+ velocity_p[i], 'kji'[i], igrid_ratio[i], v_grid_ghosts[i]))
+ s.append('{} = convert_int_rtn({});'.format(velocity_gid[i], velocity_p[i]))
+ s.append('{} = {} - {};'.format(
+ velocity_h[i], velocity_p[i], velocity_gid[i]))
+ if polynomial_interpolator:
+ s.append('{} = (int)(round({}/{}));'.format(widx[i], velocity_h[i], p_dx[i]))
else:
s.append('{} = {} + {};'.format(velocity_gid[i], 'kji'[i], velocity_ghosts))
yield ctx
@@ -544,8 +500,6 @@ 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)
local_id.declare(al,align=True,const=True)
@@ -554,19 +508,20 @@ class DirectionalAdvectionKernelGenerator(KernelCodeGenerator):
s.jumpline()
with s._align_() as al:
- compute_grid_size.declare(al,align=True)
- if has_bilevel:
- v_grid_size.declare(al,align=True)
- P_grid_ghosts.declare(al,align=True)
- V_grid_ghosts.declare(al,align=True)
+ p_grid_size.declare(al,align=True)
+ v_grid_size.declare(al,align=True)
+ p_compute_grid_size.declare(al,align=True)
+ v_compute_grid_size.declare(al,align=True)
+ p_grid_ghosts.declare(al,align=True)
+ v_grid_ghosts.declare(al,align=True)
al.jumpline()
- dx.declare(al, align=True)
- inv_dx.declare(al,align=True)
- if has_bilevel:
- v_dx.declare(al, align=True)
- v_inv_dx.declare(al,align=True)
- s.jumpline()
- xmin.declare(al,align=True)
+ p_dx.declare(al, align=True)
+ v_dx.declare(al, align=True)
+ p_inv_dx.declare(al,align=True)
+ v_inv_dx.declare(al,align=True)
+ al.jumpline()
+ p_xmin.declare(al,align=True)
+ v_xmin.declare(al,align=True)
s.jumpline()
with s._align_() as al:
@@ -576,7 +531,10 @@ class DirectionalAdvectionKernelGenerator(KernelCodeGenerator):
if is_cached:
with s._align_() as al:
- advec_ghosts.declare(al, const=True, align=True)
+ if is_bilevel:
+ grid_ratio.declare(al, const=True, align=True)
+ igrid_ratio.declare(al, const=True, align=True)
+ V_cache_ghosts.declare(al, const=True, align=True)
s.jumpline()
with s._align_() as al:
@@ -590,13 +548,15 @@ class DirectionalAdvectionKernelGenerator(KernelCodeGenerator):
kmax.declare(s)
s.decl_aligned_vars(position_gid, velocity_gid)
- if has_bilevel:
- s.decl_aligned_vars(velocity_gid_pos, velocity_h)
+ if is_bilevel:
+ s.decl_aligned_vars(velocity_p, velocity_h)
+ if polynomial_interpolator:
+ s.decl_aligned_vars(widx)
with contextlib.nested(*nested_loops):
s.jumpline()
- if is_cached and not has_bilevel:
+ if is_cached and not is_bilevel:
if tuning_mode:
loop = 'int {i}={Lx}; {i}<{N}; {i}+={gsize}'.format(
i='idx', N=V_cache_width,
@@ -611,89 +571,80 @@ class DirectionalAdvectionKernelGenerator(KernelCodeGenerator):
s.append(code)
code = 'wait_group_events(1, &event);'
s.append(code)
- s.jumpline()
- elif has_bilevel and not velocity_cache_full_length:
+ elif is_bilevel:
# We must load velocity cache line on the fly
s.comment("Load velocity in cache with linear interpolation")
- mesh_ratio = 1
- s_mesh_size = position_mesh_info['global_mesh']['compute_resolution'].value[0]
- v_mesh_size = velocity_mesh_info['global_mesh']['compute_resolution'].value[0]
- if (s_mesh_size%v_mesh_size==0):
- mesh_ratio = s_mesh_size / v_mesh_size
- with s._if_('k%{}==0'.format(mesh_ratio)):
- s.append('{} = convert_int_rtn(convert_{}({}+1)*{}*{});'.format(line_offset_for_v, ftype, line_offset, dx[0], v_inv_dx[0]))
- with s._for_('int {i}={fval}; {i}<{N} && ({i}+{o})<{Nv}; {i}+={gsize}'.format(
- i=velocity_ix, fval=local_id[0], gsize=local_size[0],
- N=Vc_shape[0],o=line_offset_for_v, Nv=v_grid_size[0]),
- unroll=False):
- if work_dim==3:
- is_first = True
- for ii,jj in [(ii,jj) for ii in range(2) for jj in range(2)]:
- s.append('{}[{}] {} ({}{})*({}{})*{};'.format(
- Vc, velocity_ix, " =" if is_first else "+=",
- "1.0-" if ii==0 else " ",velocity_h[2],
- "1.0-" if jj==0 else " ",velocity_h[1],
- s.vload(n=nparticles,ptr=line_velocity,offset='({})*{}+({})*{}+({}+{})*{}'.format(
- ii, velocity_strides[2],
- jj, velocity_strides[1],
- line_offset_for_v,
- velocity_ix, velocity_strides[0]))))
- is_first = False
- elif work_dim==2:
- for ii in range(2):
- s.append('{}[{}] {} ({}{})*{};'.format(
- Vc, velocity_ix, " =" if is_first else "+=",
- "1.0-" if jj==0 else " ",velocity_h[1],
- s.vload(n=nparticles,ptr=line_velocity,offset='({}+{})*{}+({}+{})*{}'.format(
- jj, velocity_strides[1],
- line_offset_for_v,
- velocity_ix, velocity_strides[0]))))
- is_first = False
- s.barrier(_local=True)
- s.jumpline()
+ loop = 'int {i}={Lx}; {i}<{N}; {i}+={gsize}'.format(
+ i='idx', N=V_cache_width,
+ Lx=local_id[0], gsize=local_size[0])
+ with s._for_(loop):
+ s.append('{} = ({}+{}-{})*{} + {};'.format(
+ velocity_p[i], line_offset, 'idx', V_cache_ghosts, igrid_ratio[i], v_grid_ghosts[i]))
+ s.append('{} = convert_int_rtn({});'.format(velocity_gid[i], velocity_p[i]))
+ s.append('{} = {} - {};'.format(
+ velocity_h[i], velocity_p[i], velocity_gid[i]))
+ is_first = True
+ if polynomial_interpolator:
+ s.append('{} = (int)(round({}/{}));'.format(widx[i], velocity_h[i], p_dx[i]))
+ for idx in np.ndindex(*polynomial_interpolator.n):
+ idx = idx[::-1]
+ offset = '+'.join('({}{:+})*{}'.format(gidi, idxi-gi, vsi)
+ for (gidi,idxi,gi,vsi) in zip(velocity_gid, idx, polynomial_interpolator.ghosts, velocity_strides))
+ value = '{}[{}]'.format(velocity, offset)
+ weight = '{}{}{}'.format(interpolation_weights,
+ ''.join('[{}]'.format(widx[ii]) for ii in range(work_dim)[::-1]),
+ ''.join('[{}]'.format(ii) for ii in idx[::-1]))
+ rhs = '{}*{}'.format(weight, value)
+ lhs = '{}[{}]'.format(Vc, 'idx')
+ code = '{} {} {};'.format(lhs, ' =' if is_first else '+=', rhs)
+ s.append(code)
+ is_first=False
+ else:
+ for idx in it.product(range(2), repeat=work_dim):
+ idx = idx[::-1]
+ offset = '+'.join('({}+{})*{}'.format(gidi, idxi, vsi)
+ for (gidi,idxi,vsi) in zip(velocity_gid, idx, velocity_strides))
+ value = '{}[{}]'.format(velocity, offset)
+ weight = '*'.join('({}{})'.format('1.0-' if (idxi==0) else ' ', vhi)
+ for (idxi, vhi) in zip(idx, velocity_h))
+ rhs = '{}*{}'.format(weight, value)
+ lhs = '{}[{}]'.format(Vc, 'idx')
+ code = '{} {} {};'.format(lhs, ' =' if is_first else '+=', rhs)
+ s.append(code)
+ is_first=False
+ s.barrier(_local=True)
+ s.jumpline()
with s._if_(is_active()):
pid.declare(s, init='{} + {}*{} + {}'.format(line_offset, local_id[0],
npart, poffset))
X.declare(s, init='convert_{}(min({},{}-1))*{}'.format(part_ftype, pid,
- compute_grid_size[0], dx[0]));
+ p_compute_grid_size[0], p_dx[0]));
s.jumpline()
- # Modify for bilevel interpolation : in this case, the line
- # offset is set to 0 because cache length equals local resolution
- if has_bilevel:
- rk_args = {
- 'X': X,
- 'dt': dt,
- 'inv_dx': v_inv_dx[0],
- 'line_offset': '0' if velocity_cache_full_length else line_offset_for_v,
- 'field_infos': '&{}'.format(field_infos),
- }
- else:
- rk_args = {
- 'X': X,
- 'dt': dt,
- 'inv_dx': inv_dx[0],
- 'line_offset': line_offset,
- 'field_infos': '&{}'.format(field_infos),
+ rk_args = {
+ 'X': X,
+ 'dt': dt,
+ 'inv_dx': p_inv_dx[0],
+ 'line_offset': line_offset,
}
if is_cached:
- rk_args['line_velocity'] = Vc
+ rk_args['line_velocity'] = Vc
else:
- rk_args['line_velocity'] = line_velocity
- rk_args['line_velocity'] = '{}+{}'.format(rk_args['line_velocity'], min_ghosts)
+ rk_args['line_velocity'] = line_velocity
+ rk_args['line_velocity'] = '{}+{}'.format(rk_args['line_velocity'], V_cache_ghosts)
if is_periodic:
- rk_args['line_width'] = compute_grid_size[0]
+ rk_args['line_width'] = p_compute_grid_size[0]
call = runge_kutta(**rk_args)
code = '{} = {};'.format(X(), call)
s.append(code)
if s.offset_by_xmin:
- s.append('{} += {};'.format(X(), xmin))
+ s.append('{} += {};'.format(X(), p_xmin))
s.jumpline()
s.vstore_if(cond=is_last,
- scalar_cond=lambda i: '{} < {}'.format(pid[i], compute_grid_size[0]),
+ scalar_cond=lambda i: '{} < {}'.format(pid[i], p_compute_grid_size[0]),
n=nparticles, ptr=line_position,
offset=local_id[0], data=X, offset_is_ftype=False,
use_short_circuit=use_short_circuit)
@@ -703,7 +654,7 @@ if __name__ == '__main__':
from hysop.backend.device.opencl import cl
from hysop.backend.device.codegen.base.test import _test_mesh_info, _test_typegen
- work_dim=3
+ work_dim=2
ghosts=(0,0,0)
vresolution=(128,64,32)
presolution=(256,128,64)
@@ -714,17 +665,22 @@ if __name__ == '__main__':
(_,vmesh_info) = _test_mesh_info('velocity_mesh_info',tg,work_dim,ghosts,vresolution)
(_,pmesh_info) = _test_mesh_info('position_mesh_info',tg,work_dim,ghosts,presolution)
+ grid_ratio = np.asarray([2.0,3.0,4.0][:work_dim], dtype=np.float64)
+ PI = PolynomialInterpolator(dim=work_dim, deg=3, fd=2, verbose=False)
+ polynomial_interpolator = PI.generate_subgrid_interpolator(grid_ratio=tuple(grid_ratio.astype(int)))
+
dak = DirectionalAdvectionKernelGenerator(typegen=tg, ftype=tg.fbtype,
work_dim=work_dim,
rk_scheme=ExplicitRungeKutta('Euler'),
vboundary=(BoundaryCondition.PERIODIC,BoundaryCondition.PERIODIC),
is_cached=True,
- min_ghosts=8,
+ velocity_cfl=1.42,
symbolic_mode=True,
relative_velocity=0.66,
- tuning_mode = True,
- #nparticles=4, # MONOLEVEL TEST
- nparticles=1, is_bilevel=(256,64,32), # BILEVEL TEST
+ tuning_mode = False,
+ nparticles=4,
+ grid_ratio=grid_ratio,
+ polynomial_interpolator = polynomial_interpolator,
known_vars=dict(
V_mesh_info=vmesh_info,
P_mesh_info=pmesh_info,
diff --git a/hysop/backend/device/opencl/autotunable_kernels/advection_dir.py b/hysop/backend/device/opencl/autotunable_kernels/advection_dir.py
index 654ec2be8..f5268ea2c 100644
--- a/hysop/backend/device/opencl/autotunable_kernels/advection_dir.py
+++ b/hysop/backend/device/opencl/autotunable_kernels/advection_dir.py
@@ -2,6 +2,7 @@ from hysop.tools.numpywrappers import npw
from hysop.tools.types import check_instance
from hysop.constants import AutotunerFlags, DirectionLabels, BoundaryCondition
from hysop.numerics.odesolvers.runge_kutta import ExplicitRungeKutta
+from hysop.numerics.interpolation.polynomial import PolynomialSubgridInterpolator
from hysop.backend.device.opencl import cl, clTools
from hysop.backend.device.codegen.kernels.directional_advection import DirectionalAdvectionKernelGenerator
@@ -13,7 +14,8 @@ class OpenClAutotunableDirectionalAdvectionKernel(OpenClAutotunableKernel):
def autotune(self, direction, time_integrator, velocity_cfl,
velocity, position, relative_velocity,
- hardcode_arrays, is_bilevel, **kwds):
+ hardcode_arrays, grid_ratio,
+ polynomial_interpolator, **kwds):
"""Autotune this kernel with specified configuration."""
precision = self.typegen.dtype
@@ -31,11 +33,18 @@ class OpenClAutotunableDirectionalAdvectionKernel(OpenClAutotunableKernel):
raise ValueError(msg)
if not isinstance(velocity_cfl, float) or velocity_cfl<=0.0:
raise ValueError('Invalid velocity_cfl value.')
+ check_instance(polynomial_interpolator, PolynomialSubgridInterpolator)
+
+ check_instance(grid_ratio, npw.ndarray, minval=1.0, allow_none=True)
+ is_bilevel = (grid_ratio>1.0).any()
+ periodicity = velocity.periodicity.astype(npw.int32)
+
+
+ min_velocity_ghosts = DirectionalAdvectionKernelGenerator.min_velocity_ghosts(dim, velocity_cfl, grid_ratio, periodicity)
+ velocity_cache_ghosts = DirectionalAdvectionKernelGenerator.advection_cache_ghosts(velocity_cfl, grid_ratio[-1])
- min_ghosts = DirectionalAdvectionKernelGenerator.min_ghosts(dim, velocity_cfl)
- cache_ghosts = min_ghosts[-1]
- self.check_cartesian_field(velocity, min_ghosts=min_ghosts, nb_components=1)
- if (is_bilevel is not None):
+ self.check_cartesian_field(velocity, min_ghosts=min_velocity_ghosts, nb_components=1)
+ if is_bilevel:
self.check_cartesian_field(position, nb_components=1,
compute_resolution=position.compute_resolution)
else:
@@ -53,10 +62,10 @@ class OpenClAutotunableDirectionalAdvectionKernel(OpenClAutotunableKernel):
assert isinstance(Vr, (str,float))
ftype = clTools.dtype_to_ctype(precision)
- name = 'directional_advection{}_{}__{}_{}_{}g__{}'.format(
+ name = 'directional_advection{}_{}__{}_{}__{}'.format(
"" if (is_bilevel is None) else "_bilevel_linear",
DirectionLabels[direction],
- time_integrator.name(), ftype, cache_ghosts, abs(hash(Vr)))
+ time_integrator.name(), ftype, abs(hash(Vr)))
vboundaries = (velocity.global_lboundaries[-1], velocity.global_rboundaries[-1])
@@ -94,13 +103,16 @@ class OpenClAutotunableDirectionalAdvectionKernel(OpenClAutotunableKernel):
position.mesh)
return super(OpenClAutotunableDirectionalAdvectionKernel, self).autotune(name=name,
- rk_scheme=time_integrator, kernel_args=kernel_args, cache_ghosts=cache_ghosts,
+ rk_scheme=time_integrator, kernel_args=kernel_args,
+ cache_ghosts=velocity_cache_ghosts, velocity_cfl=velocity_cfl,
vboundaries=vboundaries, precision=precision, ftype=ftype,
mesh_info_vars=mesh_info_vars, work_dim=dim,
work_size=position.compute_resolution[::-1],
hardcode_arrays=hardcode_arrays, known_args=known_args,
offset_dtype=offset_dtype, strides_dtype=strides_dtype,
- isolation_params=isolation_params, Vr=Vr, is_bilevel=is_bilevel,
+ isolation_params=isolation_params, Vr=Vr,
+ grid_ratio=grid_ratio, is_bilevel=is_bilevel,
+ polynomial_interpolator=polynomial_interpolator,
**kwds)
@@ -148,10 +160,7 @@ class OpenClAutotunableDirectionalAdvectionKernel(OpenClAutotunableKernel):
max_workitem_workload = [1,16,16]
nparticles_options = [1,2,4,8,16]
- if (extra_kwds['is_bilevel'] is not None):
- # Currently bilevel codegen is not vectorized
- nparticles_options = [1]
- max_workitem_workload = [1,1,1]
+ if extra_kwds['is_bilevel']:
caching_options = [True]
extra_kwds['max_work_load'] = max_workitem_workload[:work_dim]
@@ -190,7 +199,9 @@ class OpenClAutotunableDirectionalAdvectionKernel(OpenClAutotunableKernel):
cache_ghosts = extra_kwds['cache_ghosts']
known_args = extra_kwds['known_args']
Vr = extra_kwds['Vr']
- is_bilevel = extra_kwds['is_bilevel']
+ grid_ratio = extra_kwds['grid_ratio']
+ velocity_cfl = extra_kwds['velocity_cfl']
+ polynomial_interpolator = extra_kwds['polynomial_interpolator']
## Get compile time OpenCL known variables
known_vars = super(OpenClAutotunableDirectionalAdvectionKernel, self).generate_kernel_src(
@@ -213,27 +224,24 @@ class OpenClAutotunableDirectionalAdvectionKernel(OpenClAutotunableKernel):
rk_scheme=rk_scheme,
symbolic_mode=self.symbolic_mode,
ftype=ftype,
- min_ghosts=cache_ghosts,
+ velocity_cfl=velocity_cfl,
tuning_mode=tuning_mode,
- is_bilevel=extra_kwds['is_bilevel'],
+ grid_ratio=grid_ratio,
+ polynomial_interpolator=polynomial_interpolator,
known_vars=known_vars, **extra_parameters)
kernel_name = codegen.name
kernel_src = str(codegen)
- ## Check if cache would fit
if (local_work_size is not None):
self.check_cache(codegen.required_workgroup_cache_size(local_work_size)[2])
- if (is_bilevel is not None):
- self.check_cache(codegen.required_workgroup_velocity_cache_size())
return (kernel_name, kernel_src)
def hash_extra_kwds(self, extra_kwds):
"""Hash extra_kwds dictionnary for caching purposes."""
- kwds = ('rk_scheme', 'ftype', 'work_dim', 'Vr',
- 'vboundaries', 'cache_ghosts', 'work_size',
- 'known_args')
+ kwds = ('rk_scheme', 'ftype', 'work_dim', 'Vr', 'grid_ratio',
+ 'vboundaries', 'cache_ghosts', 'work_size', 'known_args')
return self.custom_hash(*tuple(extra_kwds[kwd] for kwd in kwds),
mesh_info_vars=extra_kwds['mesh_info_vars'])
diff --git a/hysop/backend/device/opencl/operator/directional/advection_dir.py b/hysop/backend/device/opencl/operator/directional/advection_dir.py
index 51cf8b3b8..0ada2c118 100644
--- a/hysop/backend/device/opencl/operator/directional/advection_dir.py
+++ b/hysop/backend/device/opencl/operator/directional/advection_dir.py
@@ -96,15 +96,12 @@ class OpenClDirectionalAdvection(DirectionalAdvectionBase, OpenClDirectionalOper
kwds['velocity'] = self.dvelocity
kwds['position'] = self.dposition
kwds['relative_velocity'] = self.relative_velocity
-
- if (self.is_bilevel is None):
- kwds['is_bilevel'] = None
- else:
- kwds['is_bilevel'] = self.is_bilevel
+ kwds['grid_ratio'] = self.discrete_grid_ratio
kwds['direction'] = self.splitting_direction
kwds['velocity_cfl'] = self.velocity_cfl
kwds['time_integrator'] = self.time_integrator
+ kwds['polynomial_interpolator'] = self.subgrid_interpolator
(advec_kernel, args_dict) = kernel.autotune(force_verbose=self._force_autotuner_verbose,
force_debug=self._force_autotuner_debug, hardcode_arrays=True, **kwds)
diff --git a/hysop/backend/host/python/operator/directional/advection_dir.py b/hysop/backend/host/python/operator/directional/advection_dir.py
index f746d9bb6..2a2584eed 100644
--- a/hysop/backend/host/python/operator/directional/advection_dir.py
+++ b/hysop/backend/host/python/operator/directional/advection_dir.py
@@ -71,7 +71,10 @@ class PythonDirectionalAdvection(DirectionalAdvectionBase, HostDirectionalOperat
if self.is_inplace:
self.dstmp = work.get_buffer(self, 'stmp', handle=True)
- assert (self.is_bilevel is None), "Python bilevel advection has not been implemented yet."
+ if self.is_bilevel:
+ msg="Python bilevel advection has not been implemented yet."
+ raise NotImplementedError(msg)
+
self._prepare_apply()
def _prepare_apply(self):
diff --git a/hysop/fields/cartesian_discrete_field.py b/hysop/fields/cartesian_discrete_field.py
index a7457025f..7330ae255 100644
--- a/hysop/fields/cartesian_discrete_field.py
+++ b/hysop/fields/cartesian_discrete_field.py
@@ -967,12 +967,11 @@ class CartesianDiscreteScalarFieldView(CartesianDiscreteScalarFieldViewContainer
def clone(self, name=None, pretty_name=None,
var_name=None, latex_name=None, tstate=None):
"""
- Create a new temporary DiscreteScalarField and allocate it
+ Create a DiscreteScalarField and allocate it
like the current object, possibly on a different backend.
This should only be used for debugging and testing purpose.
- The generated discrete field is not registered to the continuous
- field.
+ The generated discrete field is registered to a cloned continuous field.
"""
from hysop.tools.sympy_utils import subscript
default_name='{}__{}'.format(self.name, self._dfield._clone_id)
@@ -988,14 +987,19 @@ class CartesianDiscreteScalarFieldView(CartesianDiscreteScalarFieldViewContainer
latex_name = first_not_None(latex_name, name, default_lname)
name = first_not_None(name, default_name)
+ field = self._dfield._field
+ field = field.field_like(name=field.name+'__{}'.format(self._dfield._clone_id))
+
+ topology = self._dfield._topology.topology_like()
+
dfield = CartesianDiscreteScalarField(name=name,
pretty_name=pretty_name,
latex_name=latex_name,
var_name=var_name,
- field=self._dfield._field,
- topology=self._dfield._topology,
+ field=field,
+ topology=topology,
init_topology_state=tstate,
- register_discrete_field=False)
+ register_discrete_field=True)
dfield.copy(self)
return dfield
diff --git a/hysop/fields/continuous_field.py b/hysop/fields/continuous_field.py
index fec70928e..7f1ba42b6 100644
--- a/hysop/fields/continuous_field.py
+++ b/hysop/fields/continuous_field.py
@@ -74,9 +74,9 @@ class FieldContainerI(TaggedObject):
@abstractmethod
def tmp_like(self, name, **kwds):
- """Create a temporaty field like self, possibly altered."""
+ """Create a temporary field like self, possibly altered."""
pass
-
+
@abstractmethod
def fields(self):
"""Return all unique scalar fields contained in this field container."""
@@ -759,14 +759,14 @@ class ScalarField(NamedScalarContainerI, FieldContainerI):
check_instance(topology_state, TopologyState)
if (topology not in self.discrete_fields):
- field = topology.discretize(self)
- check_instance(field, DiscreteField)
+ dfield = topology.discretize(self)
+ check_instance(dfield, DiscreteField)
if not self.is_tmp:
assert (topology in self.discrete_fields)
- assert (self.discrete_fields[topology] is field)
+ assert (self.discrete_fields[topology] is dfield)
else:
- field = self.discrete_fields[topology]
- return field.view(topology_state)
+ dfield = self.discrete_fields[topology]
+ return dfield.view(topology_state)
def _get_dtype(self):
"""Return the default allocation dtype of this ScalarField."""
diff --git a/hysop/numerics/interpolation/polynomial.py b/hysop/numerics/interpolation/polynomial.py
index 70e83e44b..486c47628 100644
--- a/hysop/numerics/interpolation/polynomial.py
+++ b/hysop/numerics/interpolation/polynomial.py
@@ -62,6 +62,7 @@ class PolynomialInterpolator(object):
spi = str(pi)
if spi.startswith('LINEAR'):
deg=1
+ fd = 2
elif spi.startswith('CUBIC'):
deg=3
elif spi.startswith('QUINTIC'):
@@ -520,7 +521,7 @@ class PolynomialSubgridInterpolator(object):
Same data type as W.
"""
check_instance(interpolator, PolynomialInterpolator)
- check_instance(grid_ratio, tuple, size=interpolator.dim, minval=1)
+ check_instance(grid_ratio, tuple, values=int, size=interpolator.dim, minval=1)
p = interpolator.p
n = interpolator.n
diff --git a/hysop/operator/base/advection_dir.py b/hysop/operator/base/advection_dir.py
index 0b21c7a4c..0e67f8803 100644
--- a/hysop/operator/base/advection_dir.py
+++ b/hysop/operator/base/advection_dir.py
@@ -7,6 +7,7 @@ from hysop.fields.continuous_field import Field, ScalarField
from hysop.numerics.remesh.remesh import RemeshKernel
from hysop.numerics.odesolvers.runge_kutta import ExplicitRungeKutta, Euler, RK2, RK3, RK4
+from hysop.numerics.interpolation.polynomial import PolynomialInterpolation, PolynomialInterpolator
from hysop.topology.cartesian_descriptor import CartesianTopologyDescriptors
from hysop.core.memory.memory_request import MemoryRequest
@@ -17,13 +18,13 @@ class DirectionalAdvectionBase(object):
__default_method = {
TimeIntegrator: Euler,
- Interpolation: Interpolation.LINEAR,
+ Interpolation: PolynomialInterpolation.LINEAR, #Interpolation.LINEAR,
Remesh: Remesh.L2_1,
}
__available_methods = {
TimeIntegrator: InstanceOf(ExplicitRungeKutta),
- Interpolation: Interpolation.LINEAR,
+ Interpolation: (InstanceOf(Interpolation), InstanceOf(PolynomialInterpolation)),
Remesh: (InstanceOf(Remesh), InstanceOf(RemeshKernel)),
}
@@ -201,13 +202,31 @@ class DirectionalAdvectionBase(object):
s_dx = s_topo.mesh.space_step
scalar_cfl = velocity_cfl * (v_dx[ndir] / s_dx[ndir])
+ assert (v_dx>=s_dx).all()
+ grid_ratio = v_dx/s_dx
+ is_bilevel = (grid_ratio>1.0).any()
+
+ if is_bilevel:
+ interp = self.interp
+ if isinstance(interp, PolynomialInterpolation):
+ assert (grid_ratio == grid_ratio.astype(npw.int32)).all(), grid_ratio
+ polynomial_interpolator = PolynomialInterpolator.build_interpolator(
+ pi=interp, dim=self.velocity.dim)
+ elif isinstance(interp, Interpolation):
+ assert interp is Interpolation.LINEAR # legacy linear interpolator
+ polynomial_interpolator = None
+ else:
+ raise NotImplementedError(interp)
+ else:
+ polynomial_interpolator = None
+
advection_ghosts = self.velocity_advection_ghosts(velocity_cfl)
- min_velocity_ghosts = npw.integer_zeros(v_dx.shape)
- has_bilevel = any(v_dx != s_dx)
- if has_bilevel:
- # add one ghost in each direction for nd interpolation
- min_velocity_ghosts[...] = 1
- min_velocity_ghosts[ndir] = advection_ghosts
+
+ min_velocity_ghosts = self.velocity.periodicity * (grid_ratio > 1.0).astype(npw.int32)
+ min_velocity_ghosts += 1
+ if (polynomial_interpolator is not None):
+ min_velocity_ghosts += polynomial_interpolator.ghosts
+ min_velocity_ghosts[ndir] += advection_ghosts
v_requirements.min_ghosts = npw.maximum(v_requirements.min_ghosts, min_velocity_ghosts).astype(HYSOP_INTEGER)
scalar_advection_ghosts = self.velocity_advection_ghosts(scalar_cfl)
@@ -243,13 +262,16 @@ class DirectionalAdvectionBase(object):
self.advection_ghosts = advection_ghosts
self.scalar_advection_ghosts = scalar_advection_ghosts
self.remesh_ghosts = remesh_ghosts
+ self.grid_ratio = grid_ratio
+ self.is_bilevel = is_bilevel
+ self.polynomial_interpolator = polynomial_interpolator
return requirements
@classmethod
def velocity_advection_ghosts(cls, velocity_cfl):
assert velocity_cfl > 0.0, 'cfl <= 0.0'
- return int(1+npw.floor(velocity_cfl))
+ return int(npw.ceil(velocity_cfl))
@classmethod
def scalars_out_cache_ghosts(cls, scalar_cfl, remesh_kernel):
@@ -272,14 +294,21 @@ class DirectionalAdvectionBase(object):
dadvected_fields_out = {ofield: self.output_discrete_fields[ofield]
for ofield in self.advected_fields_out}
+ grid_ratio = dadvected_fields_in.values()[0].topology_state.transposed(self.grid_ratio)
+ assert (grid_ratio == (dvelocity.space_step / dadvected_fields_in.values()[0].space_step)).all()
+
+
self.dvelocity = dvelocity
self.dadvected_fields_in = dadvected_fields_in
self.dadvected_fields_out = dadvected_fields_out
+ self.discrete_grid_ratio = grid_ratio
- self.is_bilevel = None
- if any(self.dvelocity.compute_resolution != \
- self.dadvected_fields_out.values()[0].compute_resolution):
- self.is_bilevel = self.dvelocity.compute_resolution
+ if (self.polynomial_interpolator is not None) and self.is_bilevel:
+ subgrid_interpolator = self.polynomial_interpolator.generate_subgrid_interpolator(
+ grid_ratio=tuple(grid_ratio.astype(int)))
+ self.subgrid_interpolator = subgrid_interpolator
+ else:
+ self.subgrid_interpolator = None
@debug
def get_work_properties(self):
diff --git a/hysop/operator/tests/test_bilevel_advection.py b/hysop/operator/tests/test_bilevel_advection.py
index 9e5297c2d..d72b27ddd 100644
--- a/hysop/operator/tests/test_bilevel_advection.py
+++ b/hysop/operator/tests/test_bilevel_advection.py
@@ -40,7 +40,7 @@ class TestBilevelAdvectionOperator(object):
size_min=None, size_max=None):
assert dim > 0
- shape = npw.asarray((16,32,64))
+ shape = npw.asarray((16,32,64))[:dim]
npw.random.shuffle(shape)
shape = tuple(shape.tolist())
shape_s = tuple(2*s for s in shape)
@@ -74,20 +74,19 @@ class TestBilevelAdvectionOperator(object):
Vin=Vin, Sin=Sin, Sout=Sout, velocity_cfl=velocity_cfl)
@classmethod
- def __velocity_init(cls, data, coords, axes):
- for i,d in enumerate(data):
- if i in axes:
- d[...] = +1.0
- else:
- d[...] = 0.0
+ def __velocity_init(cls, data, coords, component, axes):
+ if component in axes:
+ data[...] = +1.0
+ else:
+ data[...] = 0.0
@classmethod
- def __scalar_init(cls, data, coords, offsets=None):
- offsets = first_not_None(offsets, (0.0,)*len(coords[0]))
- for i,(d,coord) in enumerate(zip(data, coords)):
- d[...] = 1.0/(i+1)
- for (c, o) in zip(coord, offsets):
- d[...] *= npw.cos(c+o)
+ def __scalar_init(cls, data, coords, component, offsets=None):
+ offsets = first_not_None(offsets, (0.0,)*len(coords))
+ assert len(coords)==len(offsets)
+ data[...] = 1.0/(component+1)
+ for (c, o) in zip(coords, offsets):
+ data[...] *= npw.cos(c+o)
def _test_one(self, time_integrator, remesh_kernel,
shape, shape_s, dim,
@@ -114,7 +113,7 @@ class TestBilevelAdvectionOperator(object):
implementations = DirectionalAdvection.implementations().keys()
implementations += Advection.implementations().keys()
implementations = list(set(implementations))
- assert ref_impl in implementations
+ assert (ref_impl in implementations)
implementations.remove(ref_impl)
implementations = [ref_impl] + implementations
@@ -269,9 +268,6 @@ class TestBilevelAdvectionOperator(object):
self._test(dim=3, is_inplace=True)
print
- def test_3D(self):
- self._test(dim=3, is_inplace=True)
-
if __name__ == '__main__':
import hysop
diff --git a/hysop/operator/tests/test_directional_advection.py b/hysop/operator/tests/test_directional_advection.py
index 05e81ccff..b424ab569 100644
--- a/hysop/operator/tests/test_directional_advection.py
+++ b/hysop/operator/tests/test_directional_advection.py
@@ -81,20 +81,19 @@ class TestDirectionalAdvectionOperator(object):
Vin=Vin, Sin=Sin, Sout=Sout, velocity_cfl=velocity_cfl)
@classmethod
- def __velocity_init(cls, data, coords, axes):
- for i,d in enumerate(data):
- if i in axes:
- d[...] = +1.0
- else:
- d[...] = 0.0
+ def __velocity_init(cls, data, coords, component, axes):
+ if component in axes:
+ data[...] = +1.0
+ else:
+ data[...] = 0.0
@classmethod
- def __scalar_init(cls, data, coords, offsets=None):
- offsets = first_not_None(offsets, (0.0,)*len(coords[0]))
- for i,(d,coord) in enumerate(zip(data, coords)):
- d[...] = 1.0/(i+1)
- for (c, o) in zip(coord, offsets):
- d[...] *= npw.cos(c+o)
+ def __scalar_init(cls, data, coords, component, offsets=None):
+ offsets = first_not_None(offsets, (0.0,)*len(coords))
+ assert len(coords)==len(offsets)
+ data[...] = 1.0/(component+1)
+ for (c, o) in zip(coords, offsets):
+ data[...] *= npw.cos(c+o)
def _test_one(self, time_integrator, remesh_kernel,
shape, dim,
@@ -115,7 +114,7 @@ class TestDirectionalAdvectionOperator(object):
# Use optimal timestep, ||Vi||_inf is 1 on a per-axis basis
dt = ScalarParameter('dt', initial_value=npw.nan)
- dt.value = (0.99 * velocity_cfl) / (max(shape)-1)
+ dt.value = (0.95 * velocity_cfl) / (max(shape)-1)
ref_impl = Implementation.OPENCL
implementations = DirectionalAdvection.implementations().keys()
@@ -212,6 +211,8 @@ class TestDirectionalAdvectionOperator(object):
if is_ref:
dxk = -Vi*(k+0)*dt()
+ dsref[0].sdata[...] = npw.nan
+ dsref[1].sdata[...] = npw.inf
dsref.initialize(self.__scalar_init, offsets=dxk.tolist())
d = dsout.distance(dsref, p=2)
if npw.any(d > 1e-1):
@@ -272,13 +273,6 @@ class TestDirectionalAdvectionOperator(object):
sys.stdout.flush()
raise
- # def test_advec_1D_out_of_place(self):
- # self._test(dim=1, is_inplace=False)
- # def test_advec_2D_out_of_place(self):
- # self._test(dim=2, is_inplace=False)
- # def test_advec_3D_out_of_place(self):
- # self._test(dim=3, is_inplace=False)
-
def test_advec_1D_inplace(self):
self._test(dim=1, is_inplace=True)
def test_advec_2D_inplace(self):
@@ -290,10 +284,6 @@ class TestDirectionalAdvectionOperator(object):
self.test_advec_1D_inplace()
self.test_advec_2D_inplace()
self.test_advec_3D_inplace()
-
- #self.test_advec_1D_out_of_place()
- #self.test_advec_2D_out_of_place()
- #self.test_advec_3D_out_of_place()
print
diff --git a/hysop/tools/parameters.py b/hysop/tools/parameters.py
index c00400c3b..900a0f943 100755
--- a/hysop/tools/parameters.py
+++ b/hysop/tools/parameters.py
@@ -106,7 +106,7 @@ class CartesianDiscretization(namedtuple("CartesianDiscretization",
lboundaries = npw.empty(shape=(resolution.size,), dtype=object)
lboundaries[...] = BoundaryCondition.PERIODIC
rboundaries = lboundaries.copy()
-
+
check_instance(lboundaries, npw.ndarray, dtype=object,
size=resolution.size, values=BoundaryCondition,
allow_none=True)
--
GitLab