diff --git a/grudge/execution.py b/grudge/execution.py
index 7e0a083569c589b8cc6645873f47a425955463c2..ddda8a79d70c1a15dd36219c99f0cfa66f9e76cd 100644
--- a/grudge/execution.py
+++ b/grudge/execution.py
@@ -54,7 +54,7 @@ class ExecutionMapper(mappers.Evaluator,
self.function_registry = bound_op.function_registry
self.queue = queue
- # {{{ expression mappings -------------------------------------------------
+ # {{{ expression mappings
def map_ones(self, expr):
if expr.dd.is_scalar():
@@ -100,6 +100,49 @@ class ExecutionMapper(mappers.Evaluator,
args = [self.rec(p) for p in expr.parameters]
return self.function_registry[expr.function.name](self.queue, *args)
+ # }}}
+
+ # {{{ elementwise reductions
+
+ def _map_elementwise_reduction(self, op_name, field_expr, dd):
+ @memoize_in(self, "elementwise_%s_knl" % op_name)
+ def knl():
+ knl = lp.make_kernel(
+ "{[el, idof, jdof]: 0<=el<nelements and 0<=idof, jdof<ndofs}",
+ """
+ result[el, idof] = %s(jdof, operand[el, jdof])
+ """ % op_name,
+ default_offset=lp.auto,
+ name="elementwise_%s_knl" % op_name)
+
+ return lp.tag_inames(knl, "el:g.0,idof:l.0")
+
+ field = self.rec(field_expr)
+ discr = self.discrwb.discr_from_dd(dd)
+ assert field.shape == (discr.nnodes,)
+
+ result = discr.empty(queue=self.queue, dtype=field.dtype,
+ allocator=self.bound_op.allocator)
+ for igrp, grp in enumerate(discr.groups):
+ knl()(self.queue,
+ operand=grp.view(field),
+ result=grp.view(result))
+
+ return result
+
+ def map_elementwise_sum(self, op, field_expr):
+ return self._map_elementwise_reduction("sum", field_expr, op.dd_in)
+
+ def map_elementwise_min(self, op, field_expr):
+ return self._map_elementwise_reduction("min", field_expr, op.dd_in)
+
+ def map_elementwise_max(self, op, field_expr):
+ return self._map_elementwise_reduction("max", field_expr, op.dd_in)
+
+ # }}}
+
+ # {{{ nodal reductions
+
def map_nodal_sum(self, op, field_expr):
# FIXME: Could allow array scalars
return cl.array.sum(self.rec(field_expr)).get()[()]
@@ -112,6 +155,8 @@ class ExecutionMapper(mappers.Evaluator,
# FIXME: Could allow array scalars
return cl.array.min(self.rec(field_expr)).get()[()]
+ # }}}
+
def map_if(self, expr):
bool_crit = self.rec(expr.condition)
@@ -152,6 +197,8 @@ class ExecutionMapper(mappers.Evaluator,
return out
+ # {{{ elementwise linear operators
+
def map_ref_diff_base(self, op, field_expr):
raise NotImplementedError(
"differentiation should be happening in batched form")
@@ -202,16 +249,6 @@ class ExecutionMapper(mappers.Evaluator,
return result
- def map_elementwise_max(self, op, field_expr):
- from grudge._internal import perform_elwise_max
- field = self.rec(field_expr)
-
- out = self.discr.volume_zeros(dtype=field.dtype)
- for eg in self.discr.element_groups:
- perform_elwise_max(eg.ranges, field, out)
-
- return out
-
def map_interpolation(self, op, field_expr):
conn = self.discrwb.connection_from_dds(op.dd_in, op.dd_out)
return conn(self.queue, self.rec(field_expr)).with_queue(self.queue)
@@ -226,6 +263,8 @@ class ExecutionMapper(mappers.Evaluator,
return self.discrwb.opposite_face_connection()(
self.queue, self.rec(field_expr)).with_queue(self.queue)
+ # }}}
+
# {{{ face mass operator
def map_ref_face_mass_operator(self, op, field_expr):
@@ -281,8 +320,6 @@ class ExecutionMapper(mappers.Evaluator,
return result
- # }}}
-
def map_signed_face_ones(self, expr):
face_discr = self.discrwb.discr_from_dd(expr.dd)
assert face_discr.dim == 0
diff --git a/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py
index b277c4c99e97280e4f82be63cac8bbcbb00bee8e..6e6e259e02a576be09cbd567885f1ee2fb2802d2 100644
--- a/grudge/symbolic/mappers/__init__.py
+++ b/grudge/symbolic/mappers/__init__.py
@@ -133,6 +133,10 @@ class OperatorReducerMixin(LocalOpReducerMixin, FluxOpReducerMixin):
map_interpolation = _map_op_base
+ map_elementwise_sum = _map_op_base
+ map_elementwise_min = _map_op_base
+ map_elementwise_max = _map_op_base
+
map_nodal_sum = _map_op_base
map_nodal_min = _map_op_base
map_nodal_max = _map_op_base
@@ -182,6 +186,10 @@ class IdentityMapperMixin(LocalOpReducerMixin, FluxOpReducerMixin):
map_interpolation = map_elementwise_linear
+ map_elementwise_sum = map_elementwise_linear
+ map_elementwise_min = map_elementwise_linear
+ map_elementwise_max = map_elementwise_linear
+
map_nodal_sum = map_elementwise_linear
map_nodal_min = map_elementwise_linear
map_nodal_max = map_elementwise_linear
@@ -746,6 +754,19 @@ class StringifyMapper(pymbolic.mapper.stringifier.StringifyMapper):
def _format_op_dd(self, op):
return ":%s->%s" % (self._format_dd(op.dd_in), self._format_dd(op.dd_out))
+ # {{{ elementwise ops
+
+ def map_elementwise_sum(self, expr, enclosing_prec):
+ return "ElementwiseSum" + self._format_op_dd(expr)
+
+ def map_elementwise_max(self, expr, enclosing_prec):
+ return "ElementwiseMax" + self._format_op_dd(expr)
+
+ def map_elementwise_min(self, expr, enclosing_prec):
+ return "ElementwiseMin" + self._format_op_dd(expr)
+
+ # }}}
+
# {{{ nodal ops
def map_nodal_sum(self, expr, enclosing_prec):
diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py
index c1360b1627a79ccafe3e52405a099f550012a410..2352b189ed3686165c7b6cba2cff0ebebe35a80d 100644
--- a/grudge/symbolic/operators.py
+++ b/grudge/symbolic/operators.py
@@ -162,9 +162,26 @@ class InterpolationOperator(Operator):
interp = InterpolationOperator
-class ElementwiseMaxOperator(Operator):
+# {{{ element reduction: sum, min, max
+
+class ElementwiseReductionOperator(Operator):
+ def __init__(self, dd):
+ super(ElementwiseReductionOperator, self).__init__(dd_in=dd, dd_out=dd)
+
+
+class ElementwiseSumOperator(ElementwiseReductionOperator):
+ mapper_method = intern("map_elementwise_sum")
+
+
+class ElementwiseMinOperator(ElementwiseReductionOperator):
+ mapper_method = intern("map_elementwise_min")
+
+
+class ElementwiseMaxOperator(ElementwiseReductionOperator):
mapper_method = intern("map_elementwise_max")
+# }}}
+
# {{{ nodal reduction: sum, integral, max
@@ -677,6 +694,19 @@ def norm(p, arg, dd=None):
else:
raise ValueError("unsupported value of p")
+
+def h_max(dd=None):
+ import grudge.symbolic.primitives as prim
+ if dd is None:
+ dd = prim.DD_VOLUME
+ dd = prim.as_dofdesc(dd)
+
+ return NodalMax(dd_in=dd)(
+ ElementwiseSumOperator(dd)(
+ MassOperator(dd_in=dd)(prim.Ones(dd))
+ )
+ )
+
# }}}
# vim: foldmethod=marker
diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py
index aa00e8e678be045fa5ae3938a295af5da3f9ff41..8ff71cc40a26a838ceb7ed667db9a7eaba715552 100644
--- a/grudge/symbolic/primitives.py
+++ b/grudge/symbolic/primitives.py
@@ -654,6 +654,7 @@ def mv_normal(dd, ambient_dim, dim=None):
def normal(dd, ambient_dim, dim=None, quadrature_tag=None):
return mv_normal(dd, ambient_dim, dim).as_vector()
+
# }}}