diff --git a/pytential/__init__.py b/pytential/__init__.py index 7a4e1f4ec91fed3ae83e72396a9e703271a5908b..b5c568f6932bf24c98c2bcb5ce34f66343b26c05 100644 --- a/pytential/__init__.py +++ b/pytential/__init__.py @@ -1,5 +1,4 @@ -from __future__ import division -from __future__ import absolute_import +from __future__ import division, absolute_import __copyright__ = "Copyright (C) 2013 Andreas Kloeckner" @@ -71,7 +70,7 @@ def integral(discr, queue, x): @memoize_on_first_arg -def _norm_op(discr, num_components): +def _norm_2_op(discr, num_components): from pytential import sym, bind if num_components is not None: from pymbolic.primitives import make_sym_vector @@ -84,7 +83,20 @@ def _norm_op(discr, num_components): sym.integral(discr.ambient_dim, discr.dim, integrand)) -def norm(discr, queue, x): +@memoize_on_first_arg +def _norm_inf_op(discr, num_components): + from pytential import sym, bind + if num_components is not None: + from pymbolic.primitives import make_sym_vector + v = make_sym_vector("arg", num_components) + max_arg = sym.abs(np.dot(sym.conj(v), v)) + else: + max_arg = sym.abs(sym.var("arg"))**2 + + return bind(discr, sym.NodeSum(max_arg)) + + +def norm(discr, queue, x, p=2): from pymbolic.geometric_algebra import MultiVector if isinstance(x, MultiVector): x = x.as_vector(np.object) @@ -93,9 +105,21 @@ def norm(discr, queue, x): if isinstance(x, np.ndarray): num_components, = x.shape - norm_op = _norm_op(discr, num_components) - from math import sqrt - return sqrt(norm_op(queue, integrand=x)) + if p == 2: + norm_op = _norm_2_op(discr, num_components) + from math import sqrt + return sqrt(norm_op(queue, integrand=x)) + + elif p == np.inf or p == "inf": + norm_op = _norm_inf_op(discr, num_components) + norm_res = norm_op(queue, arg=x) + if isinstance(norm_res, np.ndarray): + return max(norm_res) + else: + return norm_res + + else: + raise ValueError("unsupported norm order: %s" % p) __all__ = ["sym", "bind"] diff --git a/pytential/symbolic/execution.py b/pytential/symbolic/execution.py index 8f30600a3e712fce04ca79d99e0bace33251e674..53826668d7d9d87599659a464d3dc12c29fbaab9 100644 --- a/pytential/symbolic/execution.py +++ b/pytential/symbolic/execution.py @@ -50,7 +50,10 @@ class EvaluationMapper(EvaluationMapperBase): # {{{ map_XXX def map_node_sum(self, expr): - return cl.array.sum(self.rec(expr.operand)).get() + return cl.array.sum(self.rec(expr.operand)).get()[()] + + def map_node_max(self, expr): + return cl.array.max(self.rec(expr.operand)).get()[()] def map_ones(self, expr): discr = self.bound_expr.get_discretization(expr.where) diff --git a/pytential/symbolic/mappers.py b/pytential/symbolic/mappers.py index 47319841fa52fc7a6085438fb87384c406f6ede4..d77464923cac33910ebeb46a3b336267a017f1d2 100644 --- a/pytential/symbolic/mappers.py +++ b/pytential/symbolic/mappers.py @@ -60,6 +60,8 @@ class IdentityMapper(IdentityMapperBase): def map_node_sum(self, expr): return type(expr)(self.rec(expr.operand)) + map_node_max = map_node_sum + def map_num_reference_derivative(self, expr): return type(expr)(expr.ref_axes, self.rec(expr.operand), expr.where) @@ -100,6 +102,7 @@ class CombineMapper(CombineMapperBase): def map_node_sum(self, expr): return self.rec(expr.operand) + map_node_max = map_node_sum map_num_reference_derivative = map_node_sum def map_int_g(self, expr): @@ -154,6 +157,8 @@ class EvaluationMapper(EvaluationMapperBase): def map_node_sum(self, expr): return componentwise(type(expr), self.rec(expr.operand)) + map_node_max = map_node_sum + def map_node_coordinate_component(self, expr): return expr @@ -463,6 +468,9 @@ class StringifyMapper(BaseStringifyMapper): def map_node_sum(self, expr, enclosing_prec): return "NodeSum(%s)" % self.rec(expr.operand, PREC_NONE) + def map_node_max(self, expr, enclosing_prec): + return "NodeMax(%s)" % self.rec(expr.operand, PREC_NONE) + def map_node_coordinate_component(self, expr, enclosing_prec): return "x%d.%s" % (expr.ambient_axis, stringify_where(expr.where)) diff --git a/pytential/symbolic/primitives.py b/pytential/symbolic/primitives.py index c4128c58c8ded141c012c492575b2c770a44e6e3..9869afe0f705a6c4b6113d68495b03c2af68e8c9 100644 --- a/pytential/symbolic/primitives.py +++ b/pytential/symbolic/primitives.py @@ -421,9 +421,7 @@ def mean_curvature(ambient_dim, dim=None, where=None): # {{{ operators -class NodeSum(Expression): - """Implements a global sum over all discretization nodes.""" - +class NodalOperation(Expression): def __new__(cls, operand): # If the constructor is handed a multivector object, return an # object array of the operator applied to each of the @@ -443,9 +441,19 @@ class NodeSum(Expression): def __getinitargs__(self): return (self.operand,) + +class NodeSum(NodalOperation): + """Implements a global sum over all discretization nodes.""" + mapper_method = "map_node_sum" +class NodeMax(NodalOperation): + """Implements a global maximum over all discretization nodes.""" + + mapper_method = "map_node_max" + + def integral(ambient_dim, dim, operand, where=None): """A volume integral of *operand*.""" diff --git a/test/test_layer_pot_identity.py b/test/test_layer_pot_identity.py index af073cc3561c05bc8a1329d638f0877fc615b68e..ded9111777bc3d84fd9204048e4ac5cbd23d8c9d 100644 --- a/test/test_layer_pot_identity.py +++ b/test/test_layer_pot_identity.py @@ -325,10 +325,10 @@ def test_identity_convergence(ctx_getter, case, visualize=False): pt.plot(error) pt.show() - l2_error_norm = norm(density_discr, queue, error) - print(key, l2_error_norm) + linf_error_norm = norm(density_discr, queue, error, p=np.inf) + print("--->", key, linf_error_norm) - eoc_rec.add_data_point(qbx.h_max, l2_error_norm) + eoc_rec.add_data_point(qbx.h_max, linf_error_norm) if visualize: from meshmode.discretization.visualization import make_visualizer