From 164304299ee37caa9b83ad897e7f6f3c0c004aa3 Mon Sep 17 00:00:00 2001 From: Mike Campbell Date: Mon, 5 Feb 2024 16:09:50 -0600 Subject: [PATCH] Add dt_geometric_factors util, and enable dt estimate for Quads/Hexs, extend tests for it. --- grudge/dt_utils.py | 34 ++++++++++++++------ test/mesh_data.py | 2 ++ test/test_dt_utils.py | 73 ++++++++++++++++++++++++++++++++++++------- 3 files changed, 88 insertions(+), 21 deletions(-) diff --git a/grudge/dt_utils.py b/grudge/dt_utils.py index 54a9343e8..0a9b223cf 100644 --- a/grudge/dt_utils.py +++ b/grudge/dt_utils.py @@ -229,20 +229,35 @@ def h_min_from_volume( def dt_geometric_factors( dcoll: DiscretizationCollection, dd: Optional[DOFDesc] = None) -> DOFArray: r"""Computes a geometric scaling factor for each cell following - [Hesthaven_2008]_, section 6.4, defined as the inradius (radius of an - inscribed circle/sphere). + [Hesthaven_2008]_, section 6.4, For simplicial elemenents, this factor is + defined as the inradius (radius of an inscribed circle/sphere). For + non-simplicial elements, a mean length measure is returned. - Specifically, the inradius for each element is computed using the following - formula from [Shewchuk_2002]_, Table 1, for simplicial cells - (triangles/tetrahedra): + Specifically, the inradius for each simplicial element is computed using the + following formula from [Shewchuk_2002]_, Table 1 (triangles, tetrahedra): .. math:: - r_D = \frac{d V}{\sum_{i=1}^{N_{faces}} F_i}, + r_D = \frac{d~V}{\sum_{i=1}^{N_{faces}} F_i}, where :math:`d` is the topological dimension, :math:`V` is the cell volume, and :math:`F_i` are the areas of each face of the cell. + For non-simplicial elements, we use the following formula for a mean + cell size measure: + + .. math:: + + r_D = \frac{2~d~V}{\sum_{i=1}^{N_{faces}} F_i}, + + where :math:`d` is the topological dimension, :math:`V` is the cell volume, + and :math:`F_i` are the areas of each face of the cell. Other valid choices + here include the shortest, longest, average of the cell diagonals, or edges. + The value returned by this routine (i.e. the cell volume divided by the + average cell face area) is bounded by the extrema of the cell edge lengths, + is straightforward to calculate regardless of element shape, and jibes well + with the foregoing calculation for simplicial elements. + :arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one. Defaults to the base volume discretization if not provided. :returns: a frozen :class:`~meshmode.dof_array.DOFArray` containing the @@ -256,11 +271,10 @@ def dt_geometric_factors( actx = dcoll._setup_actx volm_discr = dcoll.discr_from_dd(dd) + r_fac = dcoll.dim if any(not isinstance(grp, SimplexElementGroupBase) for grp in volm_discr.groups): - raise NotImplementedError( - "Geometric factors are only implemented for simplex element groups" - ) + r_fac = 2.0*r_fac if volm_discr.dim != volm_discr.ambient_dim: from warnings import warn @@ -342,7 +356,7 @@ def dt_geometric_factors( "e,ei->ei", 1/sae_i, actx.tag_axis(1, DiscretizationDOFAxisTag(), cv_i), - tagged=(FirstAxisIsElementsTag(),)) * dcoll.dim + tagged=(FirstAxisIsElementsTag(),)) * r_fac for cv_i, sae_i in zip(cell_vols, surface_areas))))) # }}} diff --git a/test/mesh_data.py b/test/mesh_data.py index 0ccc369a1..08799d09e 100644 --- a/test/mesh_data.py +++ b/test/mesh_data.py @@ -86,6 +86,7 @@ def get_mesh(self, resolution, mesh_order): class BoxMeshBuilder(MeshBuilder): ambient_dim = 2 + group_cls = None mesh_order = 1 resolutions = [4, 8, 16] @@ -100,6 +101,7 @@ def get_mesh(self, resolution, mesh_order): return mgen.generate_regular_rect_mesh( a=self.a, b=self.b, nelements_per_axis=resolution, + group_cls=self.group_cls, order=mesh_order) diff --git a/test/test_dt_utils.py b/test/test_dt_utils.py index 9322f7d28..884aa4f78 100644 --- a/test/test_dt_utils.py +++ b/test/test_dt_utils.py @@ -47,31 +47,52 @@ @pytest.mark.parametrize("name", ["interval", "box2d", "box3d"]) -def test_geometric_factors_regular_refinement(actx_factory, name): +@pytest.mark.parametrize("tpe", [False, True]) +def test_geometric_factors_regular_refinement(actx_factory, name, tpe): from grudge.dt_utils import dt_geometric_factors + import pyopencl as cl + from grudge.array_context import TensorProductArrayContext - actx = actx_factory() + if tpe: + ctx = cl.create_some_context() + queue = cl.CommandQueue(ctx) + actx = TensorProductArrayContext(queue) + else: + actx = actx_factory() # {{{ cases + from meshmode.mesh import TensorProductElementGroup + group_cls = TensorProductElementGroup if tpe else None + if name == "interval": from mesh_data import BoxMeshBuilder - builder = BoxMeshBuilder(ambient_dim=1) + builder = BoxMeshBuilder(ambient_dim=1, group_cls=group_cls) elif name == "box2d": from mesh_data import BoxMeshBuilder - builder = BoxMeshBuilder(ambient_dim=2) + builder = BoxMeshBuilder(ambient_dim=2, group_cls=group_cls) elif name == "box3d": from mesh_data import BoxMeshBuilder - builder = BoxMeshBuilder(ambient_dim=3) + builder = BoxMeshBuilder(ambient_dim=3, group_cls=group_cls) else: raise ValueError("unknown geometry name: %s" % name) # }}} + from meshmode.discretization.poly_element import \ + LegendreGaussLobattoTensorProductGroupFactory as Lgl + + from grudge.dof_desc import DISCR_TAG_BASE + dtag_to_grp_fac = { + DISCR_TAG_BASE: Lgl(builder.order) + } if tpe else None + order = None if tpe else builder.order + min_factors = [] for resolution in builder.resolutions: mesh = builder.get_mesh(resolution, builder.mesh_order) - dcoll = DiscretizationCollection(actx, mesh, order=builder.order) + dcoll = DiscretizationCollection(actx, mesh, order=order, + discr_tag_to_group_factory=dtag_to_grp_fac) min_factors.append( actx.to_numpy( op.nodal_min(dcoll, "vol", actx.thaw(dt_geometric_factors(dcoll)))) @@ -85,7 +106,8 @@ def test_geometric_factors_regular_refinement(actx_factory, name): # Make sure it works with empty meshes mesh = builder.get_mesh(0, builder.mesh_order) - dcoll = DiscretizationCollection(actx, mesh, order=builder.order) + dcoll = DiscretizationCollection(actx, mesh, order=order, + discr_tag_to_group_factory=dtag_to_grp_fac) factors = actx.thaw(dt_geometric_factors(dcoll)) # noqa: F841 @@ -151,8 +173,23 @@ def rhs(x): @pytest.mark.parametrize("dim", [1, 2]) @pytest.mark.parametrize("degree", [2, 4]) -def test_wave_dt_estimate(actx_factory, dim, degree, visualize=False): - actx = actx_factory() +@pytest.mark.parametrize("tpe", [False, True]) +def test_wave_dt_estimate(actx_factory, dim, degree, tpe, visualize=False): + + import pyopencl as cl + from grudge.array_context import TensorProductArrayContext + + if tpe: + ctx = cl.create_some_context() + queue = cl.CommandQueue(ctx) + actx = TensorProductArrayContext(queue) + else: + actx = actx_factory() + + # {{{ cases + + from meshmode.mesh import TensorProductElementGroup + group_cls = TensorProductElementGroup if tpe else None import meshmode.mesh.generation as mgen @@ -160,10 +197,24 @@ def test_wave_dt_estimate(actx_factory, dim, degree, visualize=False): b = [1, 1, 1] mesh = mgen.generate_regular_rect_mesh( a=a[:dim], b=b[:dim], - nelements_per_axis=(3,)*dim) + nelements_per_axis=(3,)*dim, + group_cls=group_cls) + assert mesh.dim == dim - dcoll = DiscretizationCollection(actx, mesh, order=degree) + from meshmode.discretization.poly_element import \ + LegendreGaussLobattoTensorProductGroupFactory as Lgl + + from grudge.dof_desc import DISCR_TAG_BASE + order = degree + dtag_to_grp_fac = None + if tpe: + order = None + dtag_to_grp_fac = { + DISCR_TAG_BASE: Lgl(degree) + } + dcoll = DiscretizationCollection(actx, mesh, order=order, + discr_tag_to_group_factory=dtag_to_grp_fac) from grudge.models.wave import WeakWaveOperator wave_op = WeakWaveOperator(dcoll, c=1)