[pyfr] 14/88: Add support for broadcasting to the DSL.
Ghislain Vaillant
ghisvail-guest at moszumanska.debian.org
Wed Nov 16 12:05:25 UTC 2016
This is an automated email from the git hooks/post-receive script.
ghisvail-guest pushed a commit to branch master
in repository pyfr.
commit 0b140d07b9ca281e57e288572aad73f9ddf9c5e2
Author: Freddie Witherden <freddie at witherden.org>
Date: Tue Apr 19 11:45:04 2016 -0700
Add support for broadcasting to the DSL.
This enables a single per-element scalar to be broadcast to all
of the points inside of an element. The support for broadcasting
is used to greatly simplify the implementation of artificial
viscosity while also reducing memory requirements.
---
pyfr/backends/base/generator.py | 30 ++++++++++---
pyfr/backends/cuda/generator.py | 17 ++++----
pyfr/backends/mic/generator.py | 15 ++++---
pyfr/backends/opencl/generator.py | 17 ++++----
pyfr/backends/openmp/generator.py | 21 ++++-----
pyfr/solvers/baseadvecdiff/elements.py | 50 ++++++----------------
.../solvers/baseadvecdiff/kernels/shocksensor.mako | 16 ++-----
pyfr/solvers/navstokes/elements.py | 4 +-
pyfr/solvers/navstokes/kernels/tflux.mako | 2 +-
9 files changed, 82 insertions(+), 90 deletions(-)
diff --git a/pyfr/backends/base/generator.py b/pyfr/backends/base/generator.py
index 0b24937..1dd8ef3 100644
--- a/pyfr/backends/base/generator.py
+++ b/pyfr/backends/base/generator.py
@@ -19,10 +19,10 @@ class Arg(object):
def __init__(self, name, spec, body):
self.name = name
- specptn = (r'(?:(in|inout|out)\s+)?' # Intent
- r'(?:(mpi|scalar|view)\s+)?' # Attrs
- r'([A-Za-z_]\w*)' # Data type
- r'((?:\[\d+\]){0,2})$') # Constant array dimensions
+ specptn = (r'(?:(in|inout|out)\s+)?' # Intent
+ r'(?:(broadcast|mpi|scalar|view)\s+)?' # Attrs
+ r'([A-Za-z_]\w*)' # Data type
+ r'((?:\[\d+\]){0,2})$') # Dimensions
dimsptn = r'(?<=\[)\d+(?=\])'
usedptn = r'(?:[^A-Za-z]|^){0}[^A-Za-z0-9]'.format(name)
@@ -43,13 +43,18 @@ class Arg(object):
self.ncdim = len(self.cdims)
# Attributes
+ self.isbroadcast = 'broadcast' in self.attrs
self.ismpi = 'mpi' in self.attrs
self.isused = bool(re.search(usedptn, body))
self.isview = 'view' in self.attrs
self.isscalar = 'scalar' in self.attrs
self.isvector = 'scalar' not in self.attrs
- # Currently scalar arguments must be of type fpdtype_t
+ # Validation
+ if self.isbroadcast and self.intent != 'in':
+ raise ValueError('Broadcast arguments must be of intent in')
+ if self.isbroadcast and self.ncdim != 0:
+ raise ValueError('Broadcast arguments can not have dimensions')
if self.isscalar and self.dtype != 'fpdtype_t':
raise ValueError('Scalar arguments must be of type fpdtype_t')
@@ -71,9 +76,15 @@ class BaseKernelGenerator(object, metaclass=ABCMeta):
self.scalargs = [v for v in sargs if v.isscalar]
self.vectargs = [v for v in sargs if v.isvector]
+ # If we are 1D ensure that none of our arguments are broadcasts
+ if ndim == 1 and any(v.isbroadcast for v in self.vectargs):
+ raise ValueError('Broadcast arguments are not supported in 1D '
+ 'kernels')
+
# If we are 2D ensure none of our arguments are views
if ndim == 2 and any(v.isview for v in self.vectargs):
- raise ValueError('View arguments are not supported for 2D kernels')
+ raise ValueError('View arguments are not supported for 2D '
+ 'kernels')
# Similarly, check for MPI matrices
if ndim == 2 and any(v.ismpi for v in self.vectargs):
@@ -101,6 +112,9 @@ class BaseKernelGenerator(object, metaclass=ABCMeta):
# Non-stacked vector or MPI type
elif self.ndim == 1 and (va.ncdim == 0 or va.ismpi):
argt.append([np.intp])
+ # Broadcast vector
+ elif va.isbroadcast:
+ argt.append([np.intp])
# Stacked vector/matrix/stacked matrix
else:
argt.append([np.intp, np.int32])
@@ -108,6 +122,10 @@ class BaseKernelGenerator(object, metaclass=ABCMeta):
# Return
return self.ndim, argn, argt
+ def needs_lsdim(self, arg):
+ return ((self.ndim == 2 and not arg.isbroadcast) or
+ (arg.ncdim > 0 and not arg.ismpi))
+
@abstractmethod
def render(self):
pass
diff --git a/pyfr/backends/cuda/generator.py b/pyfr/backends/cuda/generator.py
index 1204664..83c229d 100644
--- a/pyfr/backends/cuda/generator.py
+++ b/pyfr/backends/cuda/generator.py
@@ -66,9 +66,7 @@ class CUDAKernelGenerator(BaseKernelGenerator):
kargs.append('{0} {1.dtype}* __restrict__ {1.name}_v'
.format(const, va).strip())
- # If we are a matrix (ndim = 2) or a non-MPI stacked
- # vector then a leading (sub) dimension is required
- if self.ndim == 2 or (va.ncdim > 0 and not va.ismpi):
+ if self.needs_lsdim(va):
kargs.append('int lsd{0.name}'.format(va))
return '__global__ void {0}({1})'.format(self.name, ', '.join(kargs))
@@ -97,13 +95,16 @@ class CUDAKernelGenerator(BaseKernelGenerator):
return '{0}_v[{1}]'.format(arg.name, ix)
def _deref_arg_array_2d(self, arg):
- # Matrix name_v[lsdim*_y + _x]
- if arg.ncdim == 0:
- ix = 'lsd{}*_y + _x'.format(arg.name)
- # Stacked matrix; name_v[(_y*nv + \1)*lsdim + _x]
+ # Broadcast vector: name_v[_x]
+ if arg.isbroadcast:
+ ix = '_x'
+ # Matrix: name_v[lsdim*_y + _x]
+ elif arg.ncdim == 0:
+ ix = 'lsd{0}*_y + _x'.format(arg.name)
+ # Stacked matrix: name_v[(_y*nv + \1)*lsdim + _x]
elif arg.ncdim == 1:
ix = r'(_y*{0} + \1)*lsd{1} + _x'.format(arg.cdims[0], arg.name)
- # Doubly stacked matrix; name_v[((\1*_ny + _y)*nv + \2)*lsdim + _x]
+ # Doubly stacked matrix: name_v[((\1*_ny + _y)*nv + \2)*lsdim + _x]
else:
ix = (r'((\1*_ny + _y)*{0} + \2)*lsd{1} + _x'
.format(arg.cdims[1], arg.name))
diff --git a/pyfr/backends/mic/generator.py b/pyfr/backends/mic/generator.py
index 4608acb..5dc09de 100644
--- a/pyfr/backends/mic/generator.py
+++ b/pyfr/backends/mic/generator.py
@@ -140,9 +140,7 @@ class MICKernelGenerator(BaseKernelGenerator):
kpack.append('{0} {1.dtype}* {1.name}_v = *arg{{0}};'
.format(const, va).strip())
- # If we are a matrix (ndim = 2) or a non-MPI stacked
- # vector then a leading (sub) dimension is required
- if self.ndim == 2 or (va.ncdim > 0 and not va.ismpi):
+ if self.needs_lsdim(va):
kspec.append('long *arg{0}')
kpack.append('int lsd{0.name} = *arg{{0}};'.format(va))
@@ -200,15 +198,18 @@ class MICKernelGenerator(BaseKernelGenerator):
def _offset_arg_array_2d(self, arg):
stmts = []
- # Matrix; name + _y*lsdim + cb
- if arg.ncdim == 0:
+ # Broadcast vector: name + cb
+ if arg.isbroadcast:
+ stmts.append('{0}_v + cb'.format(arg.name))
+ # Matrix: name + _y*lsdim + cb
+ elif arg.ncdim == 0:
stmts.append('{0}_v + _y*lsd{0} + cb'.format(arg.name))
- # Stacked matrix; name + (_y*nv + <0>)*lsdim + cb
+ # Stacked matrix: name + (_y*nv + <0>)*lsdim + cb
elif arg.ncdim == 1:
stmts.extend('{0}_v + (_y*{1} + {2})*lsd{0} + cb'
.format(arg.name, arg.cdims[0], i)
for i in range(arg.cdims[0]))
- # Doubly stacked matrix; name + ((<0>*_ny + _y)*nv + <1>)*lsdim + cb
+ # Doubly stacked matrix: name + ((<0>*_ny + _y)*nv + <1>)*lsdim + cb
else:
stmts.extend('{0}_v + (({1}*_ny + _y)*{2} + {3})*lsd{0} + cb'
.format(arg.name, i, arg.cdims[1], j)
diff --git a/pyfr/backends/opencl/generator.py b/pyfr/backends/opencl/generator.py
index c1d25ee..398db7b 100644
--- a/pyfr/backends/opencl/generator.py
+++ b/pyfr/backends/opencl/generator.py
@@ -64,9 +64,7 @@ class OpenCLKernelGenerator(BaseKernelGenerator):
else:
ka.append('__global {0.dtype}* restrict {0.name}_v')
- # If we are a matrix (ndim = 2) or a non-MPI stacked
- # vector then a leading (sub) dimension is required
- if self.ndim == 2 or (va.ncdim > 0 and not va.ismpi):
+ if self.needs_lsdim(va):
ka.append('int lsd{0.name}')
# Format
@@ -98,13 +96,16 @@ class OpenCLKernelGenerator(BaseKernelGenerator):
return '{0}_v[{1}]'.format(arg.name, ix)
def _deref_arg_array_2d(self, arg):
- # Matrix name_v[lsdim*_y + _x]
- if arg.ncdim == 0:
- ix = 'lsd{}*_y + _x'.format(arg.name)
- # Stacked matrix; name_v[(_y*nv + \1)*lsdim + _x]
+ # Broadcast vector: name_v[_x]
+ if arg.isbroadcast:
+ ix = '_x'
+ # Matrix: name_v[lsdim*_y + _x]
+ elif arg.ncdim == 0:
+ ix = 'lsd{0}*_y + _x'.format(arg.name)
+ # Stacked matrix: name_v[(_y*nv + \1)*lsdim + _x]
elif arg.ncdim == 1:
ix = r'(_y*{0} + \1)*lsd{1} + _x'.format(arg.cdims[0], arg.name)
- # Doubly stacked matrix; name_v[((\1*_ny + _y)*nv + \2)*lsdim + _x]
+ # Doubly stacked matrix: name_v[((\1*_ny + _y)*nv + \2)*lsdim + _x]
else:
ix = (r'((\1*_ny + _y)*{0} + \2)*lsd{1} + _x'
.format(arg.cdims[1], arg.name))
diff --git a/pyfr/backends/openmp/generator.py b/pyfr/backends/openmp/generator.py
index 71b55df..d9d026d 100644
--- a/pyfr/backends/openmp/generator.py
+++ b/pyfr/backends/openmp/generator.py
@@ -127,9 +127,7 @@ class OpenMPKernelGenerator(BaseKernelGenerator):
kargs.append('{0} {1.dtype}* __restrict__ {1.name}_v'
.format(const, va).strip())
- # If we are a matrix (ndim = 2) or a non-MPI stacked
- # vector then a leading (sub) dimension is required
- if self.ndim == 2 or (va.ncdim > 0 and not va.ismpi):
+ if self.needs_lsdim(va):
kargs.append('int lsd{0.name}'.format(va))
return 'void {0}({1})'.format(self.name, ', '.join(kargs))
@@ -170,13 +168,13 @@ class OpenMPKernelGenerator(BaseKernelGenerator):
# Leading (sub) dimension
lsdim = 'lsd' + arg.name if not arg.ismpi else '_nx'
- # Vector name_v[_x]
+ # Vector: name_v[_x]
if arg.ncdim == 0:
ix = '_x'
- # Stacked vector; name_v[lsdim*\1 + _x]
+ # Stacked vector: name_v[lsdim*\1 + _x]
elif arg.ncdim == 1:
ix = r'{0}*\1 + _x'.format(lsdim)
- # Doubly stacked vector; name_v[lsdim*nv*\1 + lsdim*\2 + _x]
+ # Doubly stacked vector: name_v[lsdim*nv*\1 + lsdim*\2 + _x]
else:
ix = r'{0}*{1}*\1 + {0}*\2 + _x'.format(lsdim, arg.cdims[1])
@@ -185,15 +183,18 @@ class OpenMPKernelGenerator(BaseKernelGenerator):
def _offset_arg_array_2d(self, arg):
stmts = []
- # Matrix; name + _y*lsdim + cb
- if arg.ncdim == 0:
+ # Broadcast vector: name + cb
+ if arg.isbroadcast:
+ stmts.append('{0}_v + cb'.format(arg.name))
+ # Matrix: name + _y*lsdim + cb
+ elif arg.ncdim == 0:
stmts.append('{0}_v + _y*lsd{0} + cb'.format(arg.name))
- # Stacked matrix; name + (_y*nv + <0>)*lsdim + cb
+ # Stacked matrix: name + (_y*nv + <0>)*lsdim + cb
elif arg.ncdim == 1:
stmts.extend('{0}_v + (_y*{1} + {2})*lsd{0} + cb'
.format(arg.name, arg.cdims[0], i)
for i in range(arg.cdims[0]))
- # Doubly stacked matrix; name + ((<0>*_ny + _y)*nv + <1>)*lsdim + cb
+ # Doubly stacked matrix: name + ((<0>*_ny + _y)*nv + <1>)*lsdim + cb
else:
stmts.extend('{0}_v + (({1}*_ny + _y)*{2} + {3})*lsd{0} + cb'
.format(arg.name, i, arg.cdims[1], j)
diff --git a/pyfr/solvers/baseadvecdiff/elements.py b/pyfr/solvers/baseadvecdiff/elements.py
index 437a37b..11c881d 100644
--- a/pyfr/solvers/baseadvecdiff/elements.py
+++ b/pyfr/solvers/baseadvecdiff/elements.py
@@ -48,53 +48,33 @@ class BaseAdvectionDiffusionElements(BaseAdvectionElements):
shockvar=shockvar
)
- # Allocate required scratch space for artificial viscosity
- if 'flux' in self.antialias:
- self.avis_qpts = backend.matrix((self.nqpts, 1, neles),
- extent='avis_qpts', tags=tags)
- self.avis_upts = backend.matrix((nupts, 1, neles),
- aliases=self.avis_qpts, tags=tags)
- else:
- self.avis_upts = backend.matrix((nupts, 1, neles),
- extent='avis_upts', tags=tags)
-
- if nfpts >= nupts:
- self._avis_fpts = backend.matrix((nfpts, 1, neles),
- extent='avis_fpts', tags=tags)
- avis_upts_temp = backend.matrix(
- (nupts, 1, neles), aliases=self._avis_fpts, tags=tags
- )
- else:
- avis_upts_temp = backend.matrix((nupts, 1, neles),
- extent='avis_fpts', tags=tags)
- self._avis_fpts = backend.matrix(
- (nfpts, 1, neles), aliases=avis_upts_temp, tags=tags
- )
+ # Allocate space for the artificial viscosity
+ self.avis = backend.matrix((1, neles), extent='avis', tags=tags)
+
+ # Scratch space
+ tmp_upts = backend.matrix((2*nupts, 1, neles),
+ aliases=self._vect_upts, tags=tags)
+ svar_upts = tmp_upts.rslice(0, nupts)
+ modal_svar_upts = tmp_upts.rslice(nupts, 2*nupts)
+
# Extract the scalar variable to be used for shock sensing
self.kernels['shockvar'] = lambda: backend.kernel(
'shockvar', tplargs=tplargs, dims=[nupts, neles],
- u=self.scal_upts_inb, s=self.avis_upts
+ u=self.scal_upts_inb, s=svar_upts
)
# Obtain the modal coefficients
rcpvdm = np.linalg.inv(self.basis.ubasis.vdm.T)
rcpvdm = backend.const_matrix(rcpvdm, tags={'align'})
self.kernels['shockvar_modal'] = lambda: backend.kernel(
- 'mul', rcpvdm, self.avis_upts, out=avis_upts_temp
+ 'mul', rcpvdm, svar_upts, out=modal_svar_upts
)
- if 'flux' in self.antialias:
- ame_e = self.avis_qpts
- tplargs['nrow_amu'] = self.nqpts
- else:
- ame_e = self.avis_upts
- tplargs['nrow_amu'] = nupts
-
# Apply the sensor to estimate the required artificial viscosity
self.kernels['shocksensor'] = lambda: backend.kernel(
'shocksensor', tplargs=tplargs, dims=[neles],
- s=avis_upts_temp, amu_e=ame_e, amu_f=self._avis_fpts
+ s=modal_svar_upts, amu=self.avis
)
def set_backend(self, backend, nscal_upts):
@@ -153,12 +133,10 @@ class BaseAdvectionDiffusionElements(BaseAdvectionElements):
if shock_capturing == 'artificial-viscosity':
self._set_backend_art_visc(backend)
elif shock_capturing == 'none':
- self.avis_qpts = self.avis_upts = None
+ self.avis = None
else:
raise ValueError('Invalid shock capturing scheme')
def get_avis_fpts_for_inter(self, eidx, fidx):
nfp = self.nfacefpts[fidx]
-
- rcmap = [(fpidx, eidx) for fpidx in self._srtd_face_fpts[fidx][eidx]]
- return (self._avis_fpts.mid,)*nfp, rcmap
\ No newline at end of file
+ return (self.avis.mid,)*nfp, ((0, eidx),)*nfp
diff --git a/pyfr/solvers/baseadvecdiff/kernels/shocksensor.mako b/pyfr/solvers/baseadvecdiff/kernels/shocksensor.mako
index 676c162..a91053f 100644
--- a/pyfr/solvers/baseadvecdiff/kernels/shocksensor.mako
+++ b/pyfr/solvers/baseadvecdiff/kernels/shocksensor.mako
@@ -9,11 +9,9 @@
<%pyfr:kernel name='shocksensor' ndim='1'
s='in fpdtype_t[${str(nupts)}]'
- amu_e='out fpdtype_t[${str(nrow_amu)}]'
- amu_f='out fpdtype_t[${str(nfpts)}]'>
+ amu='out fpdtype_t'>
// Smoothness indicator
fpdtype_t totEn = 0.0, pnEn = 1e-15;
- fpdtype_t se0 = ${math.log10(c['s0'])};
% for i, deg in enumerate(ubdegs):
totEn += s[${i}]*s[${i}];
@@ -22,7 +20,8 @@
% endif
% endfor
- fpdtype_t se = log10(pnEn/totEn);
+ fpdtype_t se0 = ${math.log10(c['s0'])};
+ fpdtype_t se = log10(pnEn/totEn);
// Compute cell-wise artificial viscosity
fpdtype_t mu = (se < se0 - ${c['kappa']})
@@ -30,12 +29,5 @@
: ${0.5*c['max-amu']}*(1.0 + sin(${0.5*pi/c['kappa']}*(se - se0)));
mu = (se < se0 + ${c['kappa']}) ? mu : ${c['max-amu']};
- // Copy to all upts (or qpts) and fpts
-% for i in range(nrow_amu):
- amu_e[${i}] = mu;
-% endfor
-
-% for i in range(nfpts):
- amu_f[${i}] = mu;
-% endfor
+ amu = mu;
</%pyfr:kernel>
diff --git a/pyfr/solvers/navstokes/elements.py b/pyfr/solvers/navstokes/elements.py
index b5994cc..1d22eb6 100644
--- a/pyfr/solvers/navstokes/elements.py
+++ b/pyfr/solvers/navstokes/elements.py
@@ -25,11 +25,11 @@ class NavierStokesElements(BaseFluidElements, BaseAdvectionDiffusionElements):
self.kernels['tdisf'] = lambda: backend.kernel(
'tflux', tplargs=tplargs, dims=[self.nqpts, self.neles],
u=self._scal_qpts, smats=self.smat_at('qpts'),
- f=self._vect_qpts, amu=self.avis_qpts
+ f=self._vect_qpts, amu=self.avis
)
else:
self.kernels['tdisf'] = lambda: backend.kernel(
'tflux', tplargs=tplargs, dims=[self.nupts, self.neles],
u=self.scal_upts_inb, smats=self.smat_at('upts'),
- f=self._vect_upts, amu=self.avis_upts
+ f=self._vect_upts, amu=self.avis
)
diff --git a/pyfr/solvers/navstokes/kernels/tflux.mako b/pyfr/solvers/navstokes/kernels/tflux.mako
index 813e63d..a52f3b1 100644
--- a/pyfr/solvers/navstokes/kernels/tflux.mako
+++ b/pyfr/solvers/navstokes/kernels/tflux.mako
@@ -9,7 +9,7 @@
<%pyfr:kernel name='tflux' ndim='2'
u='in fpdtype_t[${str(nvars)}]'
smats='in fpdtype_t[${str(ndims)}][${str(ndims)}]'
- amu='in fpdtype_t'
+ amu='in broadcast fpdtype_t'
f='inout fpdtype_t[${str(ndims)}][${str(nvars)}]'>
// Compute the flux (F = Fi + Fv)
fpdtype_t ftemp[${ndims}][${nvars}];
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/pyfr.git
More information about the debian-science-commits
mailing list