From 5586a7443a5728d98c45d18e0239ade86d53150a Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Thu, 10 Aug 2017 15:47:11 -0500 Subject: [PATCH] Introduce ability to easily take inf norms, start using in Green test --- pytential/__init__.py | 38 ++++++++++++++++++++++++++------ pytential/symbolic/execution.py | 5 ++++- pytential/symbolic/mappers.py | 8 +++++++ pytential/symbolic/primitives.py | 14 +++++++++--- test/test_layer_pot_identity.py | 6 ++--- 5 files changed, 57 insertions(+), 14 deletions(-) diff --git a/pytential/__init__.py b/pytential/__init__.py index 7a4e1f4e..b5c568f6 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 8f30600a..53826668 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 47319841..d7746492 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 c4128c58..9869afe0 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 af073cc3..ded91117 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 -- GitLab