From 8c650a0c178e1adc6831bf6fbadd5ad671eeb705 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sun, 3 May 2020 21:17:55 -0500 Subject: [PATCH 1/3] add elementwise reductions --- grudge/execution.py | 63 +++++++++++++++++++++++------ grudge/symbolic/mappers/__init__.py | 21 ++++++++++ grudge/symbolic/operators.py | 32 ++++++++++++++- grudge/symbolic/primitives.py | 1 + 4 files changed, 103 insertions(+), 14 deletions(-) diff --git a/grudge/execution.py b/grudge/execution.py index 7e0a0835..ddda8a79 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%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 c1360b16..2352b189 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 aa00e8e6..8ff71cc4 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() + # }}} -- GitLab From 168574c1c9f7dfdc46c409ebbc48efa4d66f84fb Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sun, 3 May 2020 21:26:39 -0500 Subject: [PATCH 2/3] use h_max in advection convergence test --- grudge/function_registry.py | 1 + grudge/symbolic/operators.py | 7 +++++-- grudge/symbolic/primitives.py | 1 - test/test_grudge.py | 15 +++++++-------- 4 files changed, 13 insertions(+), 11 deletions(-) diff --git a/grudge/function_registry.py b/grudge/function_registry.py index 2f56eba4..4d521560 100644 --- a/grudge/function_registry.py +++ b/grudge/function_registry.py @@ -136,6 +136,7 @@ class CElementwiseBinaryFunction(Function): return func(arg0, arg1) from pymbolic.primitives import Variable + @memoize_in(self, "map_call_knl_%s" % func_name) def knl(): i = Variable("i") diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py index 2352b189..1a573058 100644 --- a/grudge/symbolic/operators.py +++ b/grudge/symbolic/operators.py @@ -695,17 +695,20 @@ def norm(p, arg, dd=None): raise ValueError("unsupported value of p") -def h_max(dd=None): +def h_max(ambient_dim, dim=None, dd=None): import grudge.symbolic.primitives as prim if dd is None: dd = prim.DD_VOLUME dd = prim.as_dofdesc(dd) + if dim is None: + dim = ambient_dim + return NodalMax(dd_in=dd)( ElementwiseSumOperator(dd)( MassOperator(dd_in=dd)(prim.Ones(dd)) ) - ) + )**(1.0/dim) # }}} diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py index 8ff71cc4..aa00e8e6 100644 --- a/grudge/symbolic/primitives.py +++ b/grudge/symbolic/primitives.py @@ -654,7 +654,6 @@ 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() - # }}} diff --git a/test/test_grudge.py b/test/test_grudge.py index 65cd9d89..e27599d9 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -232,7 +232,6 @@ def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, op_type, flux_type [np.linspace(-1.0, 1.0, mesh_par)], order=order) - h = 2.0 / mesh_par dim = 1 dt_factor = 1.0 elif mesh_name == "disk": @@ -250,7 +249,6 @@ def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, op_type, flux_type from meshmode.mesh.io import from_meshpy mesh = from_meshpy(mesh_info, order=1) - h = np.sqrt(mesh_par) dim = 2 dt_factor = 4 elif mesh_name.startswith("rect"): @@ -259,7 +257,6 @@ def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, op_type, flux_type mesh = generate_regular_rect_mesh(a=(-0.5,)*dim, b=(0.5,)*dim, n=(mesh_par,)*dim, order=4) - h = 1/mesh_par if dim == 2: dt_factor = 4 elif dim == 3: @@ -306,7 +303,8 @@ def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, op_type, flux_type else: final_time = 0.2 - dt = dt_factor * h/order**2 + h_max = bind(discr, sym.h_max(discr.ambient_dim))(queue) + dt = dt_factor * h_max/order**2 nsteps = (final_time // dt) + 1 dt = final_time/nsteps + 1e-15 @@ -335,11 +333,12 @@ def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, op_type, flux_type error_l2 = bind(discr, sym.norm(2, sym.var("u")-u_analytic(sym.nodes(dim))))( queue, t=last_t, u=last_u) - print(h, error_l2) - eoc_rec.add_data_point(h, error_l2) + print(h_max, error_l2) + eoc_rec.add_data_point(h_max, error_l2) - print(eoc_rec.pretty_print(abscissa_label="h", - error_label="L2 Error")) + print(eoc_rec.pretty_print( + abscissa_label="h", + error_label="L2 Error")) assert eoc_rec.order_estimate() > order -- GitLab From 6a9610a2c0e60f821e96681129f7c00a7badef37 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sun, 17 May 2020 21:10:08 -0500 Subject: [PATCH 3/3] add docs for reduction operators --- doc/conf.py | 12 ++++++------ grudge/execution.py | 2 +- grudge/symbolic/operators.py | 21 ++++++++++++++++++++- test/test_grudge.py | 2 +- 4 files changed, 28 insertions(+), 9 deletions(-) diff --git a/doc/conf.py b/doc/conf.py index 560efec5..6bb94873 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -319,12 +319,12 @@ texinfo_documents = [ # Example configuration for intersphinx: refer to the Python standard library. intersphinx_mapping = { - 'http://docs.python.org/': None, - 'http://docs.scipy.org/doc/numpy/': None, - 'http://documen.tician.de/pyopencl/': None, - 'http://documen.tician.de/modepy/': None, - 'http://documen.tician.de/pymbolic/': None, - 'http://documen.tician.de/meshmode/': None, + 'https://docs.python.org/': None, + 'https://docs.scipy.org/doc/numpy/': None, + 'https://documen.tician.de/pyopencl/': None, + 'https://documen.tician.de/modepy/': None, + 'https://documen.tician.de/pymbolic/': None, + 'https://documen.tician.de/meshmode/': None, #'http://documen.tician.de/loopy/': None, } autoclass_content = "both" diff --git a/grudge/execution.py b/grudge/execution.py index 5423eb0a..a95b15e7 100644 --- a/grudge/execution.py +++ b/grudge/execution.py @@ -123,7 +123,7 @@ class ExecutionMapper(mappers.Evaluator, result = discr.empty(queue=self.queue, dtype=field.dtype, allocator=self.bound_op.allocator) - for igrp, grp in enumerate(discr.groups): + for grp in discr.groups: knl()(self.queue, operand=grp.view(field), result=grp.view(result)) diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py index 5678d4fb..2722e161 100644 --- a/grudge/symbolic/operators.py +++ b/grudge/symbolic/operators.py @@ -43,6 +43,8 @@ Basic Operators Reductions ^^^^^^^^^^ +.. autoclass:: ElementwiseSumOperator +.. autoclass:: ElementwiseMinOperator .. autoclass:: ElementwiseMaxOperator .. autoclass:: NodalReductionOperator @@ -178,14 +180,26 @@ class ElementwiseReductionOperator(Operator): class ElementwiseSumOperator(ElementwiseReductionOperator): + """Returns a vector of DOFs with all entries on each element set + to the sum of DOFs on that element. + """ + mapper_method = intern("map_elementwise_sum") class ElementwiseMinOperator(ElementwiseReductionOperator): + """Returns a vector of DOFs with all entries on each element set + to the minimum of DOFs on that element. + """ + mapper_method = intern("map_elementwise_min") class ElementwiseMaxOperator(ElementwiseReductionOperator): + """Returns a vector of DOFs with all entries on each element set + to the maximum of DOFs on that element. + """ + mapper_method = intern("map_elementwise_max") # }}} @@ -702,7 +716,12 @@ def norm(p, arg, dd=None): raise ValueError("unsupported value of p") -def h_max(ambient_dim, dim=None, dd=None): +def h_max_from_volume(ambient_dim, dim=None, dd=None): + """Defines a characteristic length based on the volume of the elements. + This length may not be representative if the elements have very high + aspect ratios. + """ + import grudge.symbolic.primitives as prim if dd is None: dd = prim.DD_VOLUME diff --git a/test/test_grudge.py b/test/test_grudge.py index 3d92db15..53bd4975 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -334,7 +334,7 @@ def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, op_type, flux_type else: final_time = 0.2 - h_max = bind(discr, sym.h_max(discr.ambient_dim))(queue) + h_max = bind(discr, sym.h_max_from_volume(discr.ambient_dim))(queue) dt = dt_factor * h_max/order**2 nsteps = (final_time // dt) + 1 dt = final_time/nsteps + 1e-15 -- GitLab