Enable dt estimate for quads/hexes, extend tests for it

Co-authored-by: Andreas Kloeckner <[email protected]>

Author: MTCam

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)))))
 
 # }}}

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)
 
 

test/test_dt_utils.py

@@ -33,7 +33,7 @@
         [PytestPyOpenCLArrayContextFactory,
          PytestPytatoPyOpenCLArrayContextFactory])
 
-from grudge import DiscretizationCollection
+from grudge import make_discretization_collection
 
 import grudge.op as op
 
@@ -47,22 +47,26 @@
 
 
 @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
 
     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)
 
@@ -71,7 +75,7 @@ def test_geometric_factors_regular_refinement(actx_factory, name):
     min_factors = []
     for resolution in builder.resolutions:
         mesh = builder.get_mesh(resolution, builder.mesh_order)
-        dcoll = DiscretizationCollection(actx, mesh, order=builder.order)
+        dcoll = make_discretization_collection(actx, mesh, order=builder.order)
         min_factors.append(
             actx.to_numpy(
                 op.nodal_min(dcoll, "vol", actx.thaw(dt_geometric_factors(dcoll))))
@@ -85,7 +89,7 @@ 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 = make_discretization_collection(actx, mesh, order=builder.order)
     factors = actx.thaw(dt_geometric_factors(dcoll))  # noqa: F841
 
 
@@ -115,7 +119,7 @@ def test_non_geometric_factors(actx_factory, name):
     degrees = list(range(1, 8))
     for degree in degrees:
         mesh = builder.get_mesh(1, degree)
-        dcoll = DiscretizationCollection(actx, mesh, order=degree)
+        dcoll = make_discretization_collection(actx, mesh, order=degree)
         factors.append(min(dt_non_geometric_factors(dcoll)))
 
     # Crude estimate, factors should behave like 1/N**2
@@ -134,7 +138,7 @@ def test_build_jacobian(actx_factory):
     mesh = mgen.generate_regular_rect_mesh(a=[0], b=[1], nelements_per_axis=(3,))
     assert mesh.dim == 1
 
-    dcoll = DiscretizationCollection(actx, mesh, order=1)
+    dcoll = make_discretization_collection(actx, mesh, order=1)
 
     def rhs(x):
         return 3*x**2 + 2*x + 5
@@ -151,19 +155,27 @@ 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):
+@pytest.mark.parametrize("tpe", [False, True])
+def test_wave_dt_estimate(actx_factory, dim, degree, tpe, visualize=False):
     actx = actx_factory()
 
+    # {{{ cases
+
+    from meshmode.mesh import TensorProductElementGroup
+    group_cls = TensorProductElementGroup if tpe else None
+
     import meshmode.mesh.generation as mgen
 
     a = [0, 0, 0]
     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)
+    dcoll = make_discretization_collection(actx, mesh, order=degree)
 
     from grudge.models.wave import WeakWaveOperator
     wave_op = WeakWaveOperator(dcoll, c=1)