From 4128ea848a69ae8cc16ea6165f9406fe7d189d27 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sat, 9 May 2020 11:08:58 -0500 Subject: [PATCH 1/2] operators: fix einsum in ref mass operator --- grudge/execution.py | 13 +++++- grudge/function_registry.py | 1 + grudge/symbolic/operators.py | 3 +- test/test_grudge.py | 81 +++++++++++++++++++++++++----------- 4 files changed, 70 insertions(+), 28 deletions(-) diff --git a/grudge/execution.py b/grudge/execution.py index 7e0a0835..3ba817e1 100644 --- a/grudge/execution.py +++ b/grudge/execution.py @@ -174,6 +174,7 @@ class ExecutionMapper(mappers.Evaluator, default_offset=lp.auto, name="diff") knl = lp.split_iname(knl, "i", 16, inner_tag="l.0") + knl = lp.tag_array_axes(knl, "mat", "stride:auto,stride:auto") return lp.tag_inames(knl, dict(k="g.0")) in_discr = self.discrwb.discr_from_dd(op.dd_in) @@ -287,8 +288,18 @@ class ExecutionMapper(mappers.Evaluator, face_discr = self.discrwb.discr_from_dd(expr.dd) assert face_discr.dim == 0 + from meshmode.discretization.connection import \ + ChainedDiscretizationConnection all_faces_conn = self.discrwb.connection_from_dds("vol", expr.dd) - field = face_discr.empty(self.queue, allocator=self.bound_op.allocator) + if isinstance(all_faces_conn, ChainedDiscretizationConnection): + # NOTE: this can happen if oversampling on a face. however, since + # the face_id stays the same, we can safely just look at the first + # connection + all_faces_conn = all_faces_conn.connections[0] + + field = face_discr.empty(self.queue, + dtype=self.discrwb.real_dtype, + allocator=self.bound_op.allocator) field.fill(1) for grp in all_faces_conn.groups: 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 c1360b16..8c2f7fd0 100644 --- a/grudge/symbolic/operators.py +++ b/grudge/symbolic/operators.py @@ -452,8 +452,7 @@ class RefMassOperator(RefMassOperatorBase): vand_inv_t = np.linalg.inv(vand).T weights = in_element_group.weights - - return np.einsum('c,bz,acz->abc', weights, vand_inv_t, o_vand) + return np.einsum('j,ik,jk->ij', weights, vand_inv_t, o_vand) mapper_method = intern("map_ref_mass") diff --git a/test/test_grudge.py b/test/test_grudge.py index 65cd9d89..12b7b116 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -88,39 +88,70 @@ def test_inverse_metric(ctx_factory, dim): assert err < 1e-12, (i, j, err) -def test_1d_mass_mat_trig(ctx_factory): +@pytest.mark.parametrize("ambient_dim", [1, 2, 3]) +@pytest.mark.parametrize("quad_tag", [sym.QTAG_NONE, "OVSMP"]) +def test_1d_mass_mat_trig(ctx_factory, ambient_dim, quad_tag): """Check the integral of some trig functions on an interval using the mass - matrix + matrix. """ - cl_ctx = cl.create_some_context() queue = cl.CommandQueue(cl_ctx) - from meshmode.mesh.generation import generate_regular_rect_mesh - - mesh = generate_regular_rect_mesh(a=(-4*np.pi,), b=(9*np.pi,), - n=(17,), order=1) - - discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=8) - - x = sym.nodes(1) - f = bind(discr, sym.cos(x[0])**2)(queue) - ones = bind(discr, sym.Ones(sym.DD_VOLUME))(queue) + nelements = 17 + order = 4 - mass_op = bind(discr, sym.MassOperator()(sym.var("f"))) + a = -4.0 * np.pi + b = +9.0 * np.pi + true_integral = 13*np.pi/2 * (b - a)**(ambient_dim - 1) - num_integral_1 = np.dot(ones.get(), mass_op(queue, f=f)) - num_integral_2 = np.dot(f.get(), mass_op(queue, f=ones)) - num_integral_3 = bind(discr, sym.integral(sym.var("f")))(queue, f=f) - - true_integral = 13*np.pi/2 - err_1 = abs(num_integral_1-true_integral) - err_2 = abs(num_integral_2-true_integral) - err_3 = abs(num_integral_3-true_integral) + from meshmode.discretization.poly_element import QuadratureSimplexGroupFactory + dd_quad = sym.DOFDesc(sym.DTAG_VOLUME_ALL, quad_tag) + if quad_tag is sym.QTAG_NONE: + quad_tag_to_group_factory = {} + else: + quad_tag_to_group_factory = { + quad_tag: QuadratureSimplexGroupFactory(order=2*order) + } - assert err_1 < 1e-10 - assert err_2 < 1e-10 - assert err_3 < 1e-10 + from meshmode.mesh.generation import generate_regular_rect_mesh + mesh = generate_regular_rect_mesh( + a=(a,)*ambient_dim, b=(b,)*ambient_dim, + n=(nelements,)*ambient_dim, order=1) + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order, + quad_tag_to_group_factory=quad_tag_to_group_factory) + + def _get_variables_on(dd): + sym_f = sym.var("f", dd=dd) + sym_x = sym.nodes(ambient_dim, dd=dd) + sym_ones = sym.Ones(dd) + + return sym_f, sym_x, sym_ones + + sym_f, sym_x, sym_ones = _get_variables_on(sym.DD_VOLUME) + f_volm = bind(discr, sym.cos(sym_x[0])**2)(queue).get() + ones_volm = bind(discr, sym_ones)(queue).get() + + sym_f, sym_x, sym_ones = _get_variables_on(dd_quad) + f_quad = bind(discr, sym.cos(sym_x[0])**2)(queue) + ones_quad = bind(discr, sym_ones)(queue) + + mass_op = bind(discr, sym.MassOperator(dd_quad, sym.DD_VOLUME)(sym_f)) + + num_integral_1 = np.dot(ones_volm, mass_op(queue, f=f_quad).get()) + err_1 = abs(num_integral_1 - true_integral) + assert err_1 < 5.0e-10, err_1 + + num_integral_2 = np.dot(f_volm, mass_op(queue, f=ones_quad).get()) + err_2 = abs(num_integral_2 - true_integral) + assert err_2 < 5.0e-10, err_2 + + if quad_tag is sym.QTAG_NONE: + # NOTE: `integral` always makes a square mass matrix and + # `QuadratureSimplexGroupFactory` does not have a `mass_matrix` method. + num_integral_3 = bind(discr, + sym.integral(sym_f, dd=dd_quad))(queue, f=f_quad) + err_3 = abs(num_integral_3 - true_integral) + assert err_3 < 5.0e-10, err_3 @pytest.mark.parametrize("dim", [1, 2, 3]) -- GitLab From 4d5d27c2e12fb2a579a2a09c8cc7b89dbdb61d73 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 12 May 2020 14:12:19 -0500 Subject: [PATCH 2/2] ignore quad_tags in map_signed_face_ones --- grudge/execution.py | 14 ++++++-------- test/test_grudge.py | 2 +- 2 files changed, 7 insertions(+), 9 deletions(-) diff --git a/grudge/execution.py b/grudge/execution.py index 3ba817e1..761b1472 100644 --- a/grudge/execution.py +++ b/grudge/execution.py @@ -285,17 +285,15 @@ class ExecutionMapper(mappers.Evaluator, # }}} def map_signed_face_ones(self, expr): + assert expr.dd.is_trace() face_discr = self.discrwb.discr_from_dd(expr.dd) assert face_discr.dim == 0 - from meshmode.discretization.connection import \ - ChainedDiscretizationConnection - all_faces_conn = self.discrwb.connection_from_dds("vol", expr.dd) - if isinstance(all_faces_conn, ChainedDiscretizationConnection): - # NOTE: this can happen if oversampling on a face. however, since - # the face_id stays the same, we can safely just look at the first - # connection - all_faces_conn = all_faces_conn.connections[0] + # NOTE: ignore quadrature_tags on expr.dd, since we only care about + # the face_id here + all_faces_conn = self.discrwb.connection_from_dds( + sym.DD_VOLUME, + sym.DOFDesc(expr.dd.domain_tag)) field = face_discr.empty(self.queue, dtype=self.discrwb.real_dtype, diff --git a/test/test_grudge.py b/test/test_grudge.py index 12b7b116..a55f92ba 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -90,7 +90,7 @@ def test_inverse_metric(ctx_factory, dim): @pytest.mark.parametrize("ambient_dim", [1, 2, 3]) @pytest.mark.parametrize("quad_tag", [sym.QTAG_NONE, "OVSMP"]) -def test_1d_mass_mat_trig(ctx_factory, ambient_dim, quad_tag): +def test_mass_mat_trig(ctx_factory, ambient_dim, quad_tag): """Check the integral of some trig functions on an interval using the mass matrix. """ -- GitLab