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

grudge/dt_utils.py

@@ -234,20 +234,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{V}{2~\text{max}\left(F_i\right)},
+
+    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
+    maximum 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
@@ -261,11 +276,9 @@ def dt_geometric_factors(
     actx = dcoll._setup_actx
     volm_discr = dcoll.discr_from_dd(dd)
 
-    if any(not isinstance(grp, SimplexElementGroupBase)
-           for grp in volm_discr.groups):
-        raise NotImplementedError(
-            "Geometric factors are only implemented for simplex element groups"
-        )
+    tpe = any(not isinstance(grp, SimplexElementGroupBase)
+              for grp in volm_discr.groups)
+    r_fac = 0.5 if tpe else dcoll.dim
 
     if volm_discr.dim != volm_discr.ambient_dim:
         from warnings import warn
@@ -294,61 +307,89 @@ def dt_geometric_factors(
         )
     )
 
-    if actx.supports_nonscalar_broadcasting:
-        # Compute total surface area of an element by summing over the
-        # individual face areas
-        surface_areas = DOFArray(
-            actx,
-            data=tuple(
-                actx.einsum(
-                    "fej->e",
-                    tag_axes(actx, {
-                        0: DiscretizationFaceAxisTag(),
-                        1: DiscretizationElementAxisTag(),
-                        2: DiscretizationDOFAxisTag()
+    if tpe:
+        if actx.supports_nonscalar_broadcasting:
+            surface_areas = DOFArray(
+                actx,
+                data=tuple(
+                    actx.np.max(
+                        tag_axes(actx, {
+                            0: DiscretizationFaceAxisTag(),
+                            1: DiscretizationElementAxisTag(),
                         },
+                                 face_ae_i.reshape(
+                                     vgrp.mesh_el_group.nfaces,
+                                     vgrp.nelements)),
+                        axis=0)
+                    for vgrp, face_ae_i in zip(volm_discr.groups, face_areas)))
+        else:
+            el_data_per_group = []
+            for igrp, group in enumerate(volm_discr.mesh.groups):
+                nelements = group.nelements
+                nfaces = group.nfaces
+                el_face_data = face_areas[igrp].reshape(
+                    nfaces, nelements, face_areas[igrp].shape[1])
+                el_data_np = np.ascontiguousarray(
+                    np.max(actx.to_numpy(el_face_data), axis=0)[:, 0:1])
+                el_data = actx.from_numpy(el_data_np)
+                el_data = el_data.reshape(nelements)
+                el_data_per_group.append(el_data)
+            surface_areas = DOFArray(actx, tuple(el_data_per_group))
+    else:
+        if actx.supports_nonscalar_broadcasting:
+            # Compute total surface area of an element by summing over the
+            # individual face areas
+            surface_areas = DOFArray(
+                actx,
+                data=tuple(
+                    actx.einsum(
+                        "fej->e",
+                        tag_axes(actx, {
+                            0: DiscretizationFaceAxisTag(),
+                            1: DiscretizationElementAxisTag(),
+                            2: DiscretizationDOFAxisTag()
+                        },
+                                 face_ae_i.reshape(
+                                     vgrp.mesh_el_group.nfaces,
+                                     vgrp.nelements,
+                                     face_ae_i.shape[-1])),
+                        tagged=(FirstAxisIsElementsTag(),))
+
+                    for vgrp, face_ae_i in zip(volm_discr.groups, face_areas)))
+        else:
+            surface_areas = DOFArray(
+                actx,
+                data=tuple(
+                    # NOTE: Whenever the array context can't perform nonscalar
+                    # broadcasting, elementwise reductions
+                    # (like `elementwise_integral`) repeat the *same* scalar value of
+                    # the reduction at each degree of freedom. To get a single
+                    # value for the face area (per face),
+                    # we simply average over the nodes, which gives the desired result.
+                    actx.einsum(
+                        "fej->e",
                         face_ae_i.reshape(
                             vgrp.mesh_el_group.nfaces,
                             vgrp.nelements,
-                            face_ae_i.shape[-1])),
-                    tagged=(FirstAxisIsElementsTag(),))
+                            face_ae_i.shape[-1]
+                        ) / afgrp.nunit_dofs,
+                        tagged=(FirstAxisIsElementsTag(),))
 
-                for vgrp, face_ae_i in zip(volm_discr.groups, face_areas)))
-    else:
-        surface_areas = DOFArray(
-            actx,
-            data=tuple(
-                # NOTE: Whenever the array context can't perform nonscalar
-                # broadcasting, elementwise reductions
-                # (like `elementwise_integral`) repeat the *same* scalar value of
-                # the reduction at each degree of freedom. To get a single
-                # value for the face area (per face),
-                # we simply average over the nodes, which gives the desired result.
-                actx.einsum(
-                    "fej->e",
-                    face_ae_i.reshape(
-                        vgrp.mesh_el_group.nfaces,
-                        vgrp.nelements,
-                        face_ae_i.shape[-1]
-                    ) / afgrp.nunit_dofs,
-                    tagged=(FirstAxisIsElementsTag(),))
-
-                for vgrp, afgrp, face_ae_i in zip(volm_discr.groups,
-                                                  face_discr.groups,
-                                                  face_areas)
+                    for vgrp, afgrp, face_ae_i in zip(volm_discr.groups,
+                                                      face_discr.groups,
+                                                      face_areas)
+                )
             )
-        )
 
     return actx.freeze(
-            actx.tag(NameHint(f"dt_geometric_{dd.as_identifier()}"),
-                DOFArray(actx,
-                    data=tuple(
-                        actx.einsum(
-                            "e,ei->ei",
-                            1/sae_i,
-                            actx.tag_axis(1, DiscretizationDOFAxisTag(), cv_i),
-                            tagged=(FirstAxisIsElementsTag(),)) * dcoll.dim
-                        for cv_i, sae_i in zip(cell_vols, surface_areas)))))
+            actx.tag(
+                NameHint(f"dt_geometric_{dd.as_identifier()}"),
+                DOFArray(actx, data=tuple(
+                    actx.einsum(
+                        "e,ei->ei", 1/sae_i,
+                        actx.tag_axis(1, DiscretizationDOFAxisTag(), cv_i),
+                        tagged=(FirstAxisIsElementsTag(),)) * r_fac
+                    for cv_i, sae_i in zip(cell_vols, surface_areas)))))
 
 # }}}
 

grudge/models/wave.py

@@ -28,12 +28,13 @@
 
 import numpy as np
 
+from meshmode.discretization.poly_element import SimplexElementGroupBase
 from meshmode.mesh import BTAG_ALL, BTAG_NONE
 from pytools.obj_array import flat_obj_array
 
 import grudge.geometry as geo
 import grudge.op as op
-from grudge.dof_desc import DISCR_TAG_BASE, as_dofdesc
+from grudge.dof_desc import DD_VOLUME_ALL, DISCR_TAG_BASE, as_dofdesc
 from grudge.models import HyperbolicOperator
 
 
@@ -187,7 +188,12 @@ def max_characteristic_velocity(self, actx, t=None, fields=None):
 
     def estimate_rk4_timestep(self, actx, dcoll, **kwargs):
         # FIXME: Sketchy, empirically determined fudge factor
-        return 0.38 * super().estimate_rk4_timestep(actx,  dcoll, **kwargs)
+        volm_discr = dcoll.discr_from_dd(DD_VOLUME_ALL)
+        tpe = any(not isinstance(grp, SimplexElementGroupBase)
+                  for grp in volm_discr.groups)
+        fudge_factor = 0.233 if tpe else 0.38
+        return fudge_factor * super().estimate_rk4_timestep(
+            actx,  dcoll, **kwargs)
 
 # }}}
 

test/mesh_data.py

@@ -4,13 +4,14 @@
 import numpy as np
 
 import meshmode.mesh.generation as mgen
-from meshmode.mesh import Mesh
+from meshmode.mesh import Mesh, MeshElementGroup
 from meshmode.mesh.io import read_gmsh
 
 
 class MeshBuilder(ABC):
     resolutions: ClassVar[Sequence[Hashable]]
     ambient_dim: ClassVar[int]
+    group_cls: ClassVar[MeshElementGroup]
 
     @abstractmethod
     def get_mesh(
@@ -110,8 +111,13 @@ def get_mesh(self, resolution, mesh_order=4):
         return affine_map(mesh, A=np.diag([1.0, 1.0, self.aspect_ratio]))
 
 
-class _BoxMeshBuilderBase(MeshBuilder):
+class BoxMeshBuilder(MeshBuilder):
     resolutions: ClassVar[Sequence[Hashable]] = [4, 8, 16]
+
+    def __init__(self, ambient_dim=2, group_cls=None):
+        self.ambient_dim = ambient_dim
+        self.group_cls = group_cls
+
     mesh_order = 1
 
     a = (-0.5, -0.5, -0.5)
@@ -124,18 +130,19 @@ def get_mesh(self, resolution, mesh_order=4):
         return mgen.generate_regular_rect_mesh(
                 a=self.a, b=self.b,
                 nelements_per_axis=resolution,
+                group_cls=self.group_cls,
                 order=mesh_order)
 
 
-class BoxMeshBuilder1D(_BoxMeshBuilderBase):
+class BoxMeshBuilder1D(BoxMeshBuilder):
     ambient_dim = 1
 
 
-class BoxMeshBuilder2D(_BoxMeshBuilderBase):
+class BoxMeshBuilder2D(BoxMeshBuilder):
     ambient_dim = 2
 
 
-class BoxMeshBuilder3D(_BoxMeshBuilderBase):
+class BoxMeshBuilder3D(BoxMeshBuilder):
     ambient_dim = 2