From f29701705b374ede4b59bac11247c0517c756d54 Mon Sep 17 00:00:00 2001 From: Mike Campbell Date: Wed, 25 Sep 2024 12:40:56 -0500 Subject: [PATCH 1/5] Add test to demonstrate simplex integral fails with overintegration --- test/mesh_data.py | 10 ++++++++-- test/test_grudge.py | 30 ++++++++++++++++++++---------- 2 files changed, 28 insertions(+), 12 deletions(-) diff --git a/test/mesh_data.py b/test/mesh_data.py index 58c9b80a7..11a9f04d9 100644 --- a/test/mesh_data.py +++ b/test/mesh_data.py @@ -117,14 +117,20 @@ class _BoxMeshBuilderBase(MeshBuilder): a = (-0.5, -0.5, -0.5) b = (+0.5, +0.5, +0.5) + tpe: bool + + def __init__(self, tpe=False): + self._tpe = tpe + def get_mesh(self, resolution, mesh_order=4): if not isinstance(resolution, (list, tuple)): resolution = (resolution,) * self.ambient_dim - + from meshmode.mesh import TensorProductElementGroup + group_cls = TensorProductElementGroup if self._tpe else None return mgen.generate_regular_rect_mesh( a=self.a, b=self.b, nelements_per_axis=resolution, - order=mesh_order) + order=mesh_order, group_cls=group_cls) class BoxMeshBuilder1D(_BoxMeshBuilderBase): diff --git a/test/test_grudge.py b/test/test_grudge.py index b0849310b..83d49a2cd 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -155,26 +155,31 @@ def _spheroid_surface_area(radius, aspect_ratio): @pytest.mark.parametrize("name", [ - "2-1-ellipse", "spheroid", "box2d", "box3d" + "box2d-tpe", "box3d-tpe", "box2d", "box3d", "2-1-ellipse", "spheroid", ]) -def test_mass_surface_area(actx_factory, name): +@pytest.mark.parametrize("oi", [False, True]) +def test_mass_surface_area(actx_factory, name, oi): + from grudge.dof_desc import as_dofdesc, DISCR_TAG_BASE, DISCR_TAG_QUAD actx = actx_factory() + qtag = DISCR_TAG_QUAD if oi else DISCR_TAG_BASE + vol_dd_base = as_dofdesc(dof_desc.DTAG_VOLUME_ALL) + vol_dd_quad = vol_dd_base.with_discr_tag(qtag) # {{{ cases order = 4 - + tpe = name.endswith("-tpe") if name == "2-1-ellipse": builder = mesh_data.EllipseMeshBuilder(radius=3.1, aspect_ratio=2.0) surface_area = _ellipse_surface_area(builder.radius, builder.aspect_ratio) elif name == "spheroid": builder = mesh_data.SpheroidMeshBuilder() surface_area = _spheroid_surface_area(builder.radius, builder.aspect_ratio) - elif name == "box2d": - builder = mesh_data.BoxMeshBuilder2D() + elif name.startswith("box2d"): + builder = mesh_data.BoxMeshBuilder2D(tpe=tpe) surface_area = 1.0 - elif name == "box3d": - builder = mesh_data.BoxMeshBuilder3D() + elif name.startswith("box3d"): + builder = mesh_data.BoxMeshBuilder3D(tpe=tpe) surface_area = 1.0 else: raise ValueError(f"unknown geometry name: {name}") @@ -189,15 +194,20 @@ def test_mass_surface_area(actx_factory, name): for resolution in builder.resolutions: mesh = builder.get_mesh(resolution, order) - dcoll = make_discretization_collection(actx, mesh, order=order) - volume_discr = dcoll.discr_from_dd(dof_desc.DD_VOLUME_ALL) + dcoll = make_discretization_collection( + actx, mesh, + discr_tag_to_group_factory={ + DISCR_TAG_BASE: InterpolatoryEdgeClusteredGroupFactory(order), + DISCR_TAG_QUAD: QuadratureGroupFactory(3 * order) + }) + volume_discr = dcoll.discr_from_dd(vol_dd_quad) logger.info("ndofs: %d", volume_discr.ndofs) logger.info("nelements: %d", volume_discr.mesh.nelements) # {{{ compute surface area - dd = dof_desc.DD_VOLUME_ALL + dd = vol_dd_quad ones_volm = volume_discr.zeros(actx) + 1 approx_surface_area = actx.to_numpy(op.integral(dcoll, dd, ones_volm)) From 3502174a4ccde5e8dca28b4ec137ee8d09f6193f Mon Sep 17 00:00:00 2001 From: Mike Campbell Date: Wed, 25 Sep 2024 13:13:07 -0500 Subject: [PATCH 2/5] Fix up linting issues --- test/mesh_data.py | 4 ++-- test/test_grudge.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/test/mesh_data.py b/test/mesh_data.py index 11a9f04d9..26b222ecc 100644 --- a/test/mesh_data.py +++ b/test/mesh_data.py @@ -125,8 +125,8 @@ def __init__(self, tpe=False): def get_mesh(self, resolution, mesh_order=4): if not isinstance(resolution, (list, tuple)): resolution = (resolution,) * self.ambient_dim - from meshmode.mesh import TensorProductElementGroup - group_cls = TensorProductElementGroup if self._tpe else None + from meshmode.mesh import TensorProductElementGroup + group_cls = TensorProductElementGroup if self._tpe else None return mgen.generate_regular_rect_mesh( a=self.a, b=self.b, nelements_per_axis=resolution, diff --git a/test/test_grudge.py b/test/test_grudge.py index 83d49a2cd..b23543605 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -159,7 +159,7 @@ def _spheroid_surface_area(radius, aspect_ratio): ]) @pytest.mark.parametrize("oi", [False, True]) def test_mass_surface_area(actx_factory, name, oi): - from grudge.dof_desc import as_dofdesc, DISCR_TAG_BASE, DISCR_TAG_QUAD + from grudge.dof_desc import DISCR_TAG_BASE, DISCR_TAG_QUAD, as_dofdesc actx = actx_factory() qtag = DISCR_TAG_QUAD if oi else DISCR_TAG_BASE vol_dd_base = as_dofdesc(dof_desc.DTAG_VOLUME_ALL) From 6b1381dcfb3e75d694cfb79836ace026acfbce01 Mon Sep 17 00:00:00 2001 From: Mike Campbell Date: Fri, 27 Sep 2024 12:20:43 -0500 Subject: [PATCH 3/5] Use quadrature weights directly instead of mass_matrix in integral utils --- grudge/reductions.py | 38 ++++++++++++++++++++++++++++---------- 1 file changed, 28 insertions(+), 10 deletions(-) diff --git a/grudge/reductions.py b/grudge/reductions.py index 6dcde1316..6f6108462 100644 --- a/grudge/reductions.py +++ b/grudge/reductions.py @@ -295,14 +295,23 @@ def integral(dcoll: DiscretizationCollection, dd, vec) -> Scalar: :class:`~arraycontext.ArrayContainer` of them. :returns: a device scalar denoting the evaluated integral. """ - from grudge.op import _apply_mass_operator - dd = dof_desc.as_dofdesc(dd) + discr = dcoll.discr_from_dd(dd) - ones = dcoll.discr_from_dd(dd).zeros(vec.array_context) + 1.0 - return nodal_sum( - dcoll, dd, vec * _apply_mass_operator(dcoll, dd, dd, ones) + from grudge.geometry import area_element + actx = vec.array_context + area_elements = area_element( + actx, dcoll, dd=dd, + _use_geoderiv_connection=actx.supports_nonscalar_broadcasting) + qwts = DOFArray( + actx, + data=tuple( + actx.from_numpy( + np.tile(grp.quadrature_rule().weights, + (ae_i.shape[0], 1)))*ae_i + for grp, ae_i in zip(discr.groups, area_elements)) ) + return nodal_sum(dcoll, dd, vec*qwts) # }}} @@ -509,13 +518,22 @@ def elementwise_integral( raise TypeError("invalid number of arguments") dd = dof_desc.as_dofdesc(dd) + discr = dcoll.discr_from_dd(dd) - from grudge.op import _apply_mass_operator - - ones = dcoll.discr_from_dd(dd).zeros(vec.array_context) + 1.0 - return elementwise_sum( - dcoll, dd, vec * _apply_mass_operator(dcoll, dd, dd, ones) + from grudge.geometry import area_element + actx = vec.array_context + area_elements = area_element( + actx, dcoll, dd=dd, + _use_geoderiv_connection=actx.supports_nonscalar_broadcasting) + qwts = DOFArray( + actx, + data=tuple( + actx.from_numpy( + np.tile(grp.quadrature_rule().weights, + (ae_i.shape[0], 1)))*ae_i + for grp, ae_i in zip(discr.groups, area_elements)) ) + return elementwise_sum(dcoll, dd, vec*qwts) # }}} From 88a4eca3677718cac9f0eda886c537d6d9a03fab Mon Sep 17 00:00:00 2001 From: Mike Campbell Date: Wed, 2 Oct 2024 10:40:07 -0500 Subject: [PATCH 4/5] Add comment about test name after these changes --- test/test_grudge.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/test/test_grudge.py b/test/test_grudge.py index b23543605..44271116f 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -153,7 +153,9 @@ def _spheroid_surface_area(radius, aspect_ratio): e = np.sqrt(1.0 - (c/a)**2) return 2.0 * np.pi * radius**2 * (1 + (c/a)**2 / e * np.arctanh(e)) - +# This tests uses op.integral which has been changed to directly use quadrature +# weights instead of the mass matrix. We either need to reimplement this test +# to do what it was intended to do, or rename it appropriately. @pytest.mark.parametrize("name", [ "box2d-tpe", "box3d-tpe", "box2d", "box3d", "2-1-ellipse", "spheroid", ]) From 10cd445c184a6f0f755d1fabadc0f964833ea2af Mon Sep 17 00:00:00 2001 From: Mike Campbell Date: Wed, 2 Oct 2024 10:47:08 -0500 Subject: [PATCH 5/5] Smooth out ruffness --- test/test_grudge.py | 1 + 1 file changed, 1 insertion(+) diff --git a/test/test_grudge.py b/test/test_grudge.py index 44271116f..181557276 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -153,6 +153,7 @@ def _spheroid_surface_area(radius, aspect_ratio): e = np.sqrt(1.0 - (c/a)**2) return 2.0 * np.pi * radius**2 * (1 + (c/a)**2 / e * np.arctanh(e)) + # This tests uses op.integral which has been changed to directly use quadrature # weights instead of the mass matrix. We either need to reimplement this test # to do what it was intended to do, or rename it appropriately.