TPE for acoustic_pulse example

examples/euler/acoustic_pulse.py

@@ -37,6 +37,7 @@
     EulerOperator,
     InviscidWallBC
 )
+from meshmode.mesh import TensorProductElementGroup
 from grudge.shortcuts import rk4_step
 
 from meshmode.mesh import BTAG_ALL
@@ -112,7 +113,8 @@ def run_acoustic_pulse(actx,
                        final_time=1,
                        resolution=16,
                        overintegration=False,
-                       visualize=False):
+                       visualize=False,
+                       tpe=False):
 
     # eos-related parameters
     gamma = 1.4
@@ -124,16 +126,18 @@ def run_acoustic_pulse(actx,
     dim = 2
     box_ll = -0.5
     box_ur = 0.5
+    group_cls = TensorProductElementGroup if tpe else None
     mesh = generate_regular_rect_mesh(
         a=(box_ll,)*dim,
         b=(box_ur,)*dim,
-        nelements_per_axis=(resolution,)*dim)
+        nelements_per_axis=(resolution,)*dim,
+        group_cls=group_cls)
 
     from grudge import DiscretizationCollection
     from grudge.dof_desc import DISCR_TAG_BASE, DISCR_TAG_QUAD
-    from meshmode.discretization.poly_element import \
-        (default_simplex_group_factory,
-         QuadratureSimplexGroupFactory)
+    from meshmode.discretization.poly_element import (
+        InterpolatoryEdgeClusteredGroupFactory,
+        QuadratureGroupFactory)
 
     exp_name = f"fld-acoustic-pulse-N{order}-K{resolution}"
     if overintegration:
@@ -145,9 +149,8 @@ def run_acoustic_pulse(actx,
     dcoll = DiscretizationCollection(
         actx, mesh,
         discr_tag_to_group_factory={
-            DISCR_TAG_BASE: default_simplex_group_factory(
-                base_dim=mesh.dim, order=order),
-            DISCR_TAG_QUAD: QuadratureSimplexGroupFactory(2*order)
+            DISCR_TAG_BASE: InterpolatoryEdgeClusteredGroupFactory(order),
+            DISCR_TAG_QUAD: QuadratureGroupFactory(2*order)
         }
     )
 
@@ -212,7 +215,8 @@ def rhs(t, q):
 
 
 def main(ctx_factory, order=3, final_time=1, resolution=16,
-         overintegration=False, visualize=False, lazy=False):
+         overintegration=False, visualize=False, lazy=False,
+         tpe=False):
     cl_ctx = ctx_factory()
     queue = cl.CommandQueue(cl_ctx)
 
@@ -234,7 +238,7 @@ def main(ctx_factory, order=3, final_time=1, resolution=16,
         resolution=resolution,
         overintegration=overintegration,
         final_time=final_time,
-        visualize=visualize
+        visualize=visualize, tpe=tpe
     )
 
 
@@ -251,6 +255,8 @@ def main(ctx_factory, order=3, final_time=1, resolution=16,
                         help="write out vtk output")
     parser.add_argument("--lazy", action="store_true",
                         help="switch to a lazy computation mode")
+    parser.add_argument("--tpe", action="store_true",
+                        help="use tensor product elements")
     args = parser.parse_args()
 
     logging.basicConfig(level=logging.INFO)
@@ -260,4 +266,4 @@ def main(ctx_factory, order=3, final_time=1, resolution=16,
          resolution=args.resolution,
          overintegration=args.oi,
          visualize=args.visualize,
-         lazy=args.lazy)
+         lazy=args.lazy, tpe=args.tpe)

grudge/geometry/metrics.py

@@ -86,6 +86,11 @@
 register_multivector_as_array_container()
 
 
+def _has_geoderiv_connection(grp):
+    from modepy.shapes import Simplex
+    return grp.is_affine and issubclass(grp._modepy_shape_cls, Simplex)
+
+
 def _geometry_to_quad_if_requested(
         dcoll, inner_dd, dd, vec, _use_geoderiv_connection):
 
@@ -105,7 +110,7 @@ def to_quad(vec):
     return DOFArray(
             vec.array_context,
             tuple(
-                geoderiv_vec_i if megrp.is_affine else all_quad_vec_i
+                geoderiv_vec_i if _has_geoderiv_connection(megrp) else all_quad_vec_i
                 for megrp, geoderiv_vec_i, all_quad_vec_i in zip(
                     dcoll.discr_from_dd(inner_dd).mesh.groups,
                     dcoll._base_to_geoderiv_connection(inner_dd)(vec),

grudge/models/euler.py

@@ -322,12 +322,6 @@ def operator(self, t, q):
         def interp_to_quad(u):
             return op.project(dcoll, "vol", dq, u)
 
-        # Compute volume fluxes
-        volume_fluxes = op.weak_local_div(
-            dcoll, dq,
-            interp_to_quad(euler_volume_flux(dcoll, q, gamma=gamma))
-        )
-
         # Compute interior interface fluxes
         interface_fluxes = (
             sum(
@@ -357,6 +351,12 @@ def interp_to_quad(u):
             )
             interface_fluxes = interface_fluxes + bc_fluxes
 
+        # Compute volume fluxes
+        volume_fluxes = op.weak_local_div(
+            dcoll, dq,
+            interp_to_quad(euler_volume_flux(dcoll, q, gamma=gamma))
+        )
+
         return op.inverse_mass(
             dcoll,
             volume_fluxes - op.face_mass(dcoll, df, interface_fluxes)

test/mesh_data.py

@@ -4,6 +4,7 @@
 from meshmode.mesh.io import read_gmsh
 import numpy as np
 import meshmode.mesh.generation as mgen
+from meshmode.mesh import TensorProductElementGroup
 
 
 class MeshBuilder(ABC):
@@ -111,6 +112,7 @@ def get_mesh(self, resolution, mesh_order=4):
 class _BoxMeshBuilderBase(MeshBuilder):
     resolutions = [4, 8, 16]
     mesh_order = 1
+    group_cls = None
 
     a = (-0.5, -0.5, -0.5)
     b = (+0.5, +0.5, +0.5)
@@ -122,20 +124,32 @@ def get_mesh(self, resolution, mesh_order=4):
         return mgen.generate_regular_rect_mesh(
                 a=self.a, b=self.b,
                 nelements_per_axis=resolution,
-                order=mesh_order)
+                order=mesh_order, group_cls=self.group_cls)
 
 
 class BoxMeshBuilder1D(_BoxMeshBuilderBase):
     ambient_dim = 1
 
+    def __init__(self, tpe=False):
+        if tpe:
+            self.group_cls = TensorProductElementGroup
+
 
 class BoxMeshBuilder2D(_BoxMeshBuilderBase):
     ambient_dim = 2
 
+    def __init__(self, tpe=False):
+        if tpe:
+            self.group_cls = TensorProductElementGroup
+
 
 class BoxMeshBuilder3D(_BoxMeshBuilderBase):
     ambient_dim = 2
 
+    def __init__(self, tpe=False):
+        if tpe:
+            self.group_cls = TensorProductElementGroup
+
 
 class WarpedRectMeshBuilder(MeshBuilder):
     resolutions = [4, 6, 8]

test/test_grudge.py

@@ -61,10 +61,10 @@
 
 # {{{ mass operator trig integration
 
+@pytest.mark.parametrize("tpe", [True, False])
 @pytest.mark.parametrize("ambient_dim", [1, 2, 3])
-@pytest.mark.parametrize("discr_tag", [dof_desc.DISCR_TAG_BASE,
-                                       dof_desc.DISCR_TAG_QUAD])
-def test_mass_mat_trig(actx_factory, ambient_dim, discr_tag):
+@pytest.mark.parametrize("use_overint", [False, True])
+def test_mass_mat_trig(actx_factory, tpe, ambient_dim, use_overint):
     """Check the integral of some trig functions on an interval using the mass
     matrix.
     """
@@ -75,22 +75,26 @@ def test_mass_mat_trig(actx_factory, ambient_dim, discr_tag):
 
     a = -4.0 * np.pi
     b = +9.0 * np.pi
+
     true_integral = 13*np.pi/2 * (b - a)**(ambient_dim - 1)
+    discr_tag = dof_desc.DISCR_TAG_QUAD if use_overint else dof_desc.DISCR_TAG_BASE
 
-    from meshmode.discretization.poly_element import QuadratureSimplexGroupFactory
     dd_quad = dof_desc.DOFDesc(dof_desc.DTAG_VOLUME_ALL, discr_tag)
+    discr_order = order
     if discr_tag is dof_desc.DISCR_TAG_BASE:
         discr_tag_to_group_factory = {}
     else:
+        discr_order = None
         discr_tag_to_group_factory = {
-            discr_tag: QuadratureSimplexGroupFactory(order=2*order)
+            dof_desc.DISCR_TAG_BASE: InterpolatoryEdgeClusteredGroupFactory(order),
+            dof_desc.DISCR_TAG_QUAD: QuadratureGroupFactory(order=2*order)
         }
 
     mesh = mgen.generate_regular_rect_mesh(
             a=(a,)*ambient_dim, b=(b,)*ambient_dim,
             nelements_per_axis=(nel_1d,)*ambient_dim, order=1)
     dcoll = DiscretizationCollection(
-        actx, mesh, order=order,
+        actx, mesh, order=discr_order,
         discr_tag_to_group_factory=discr_tag_to_group_factory
     )
 
@@ -160,27 +164,32 @@ def _spheroid_surface_area(radius, aspect_ratio):
         return 2.0 * np.pi * radius**2 * (1 + (c/a)**2 / e * np.arctanh(e))
 
 
+@pytest.mark.parametrize("tpe", [True, False])
 @pytest.mark.parametrize("name", [
     "2-1-ellipse", "spheroid", "box2d", "box3d"
     ])
-def test_mass_surface_area(actx_factory, name):
+def test_mass_surface_area(actx_factory, tpe, name):
     actx = actx_factory()
 
     # {{{ cases
 
     order = 4
 
     if name == "2-1-ellipse":
+        if tpe:
+            pytest.skip()
         builder = mesh_data.EllipseMeshBuilder(radius=3.1, aspect_ratio=2.0)
         surface_area = _ellipse_surface_area(builder.radius, builder.aspect_ratio)
     elif name == "spheroid":
+        if tpe:
+            pytest.skip()
         builder = mesh_data.SpheroidMeshBuilder()
         surface_area = _spheroid_surface_area(builder.radius, builder.aspect_ratio)
     elif name == "box2d":
-        builder = mesh_data.BoxMeshBuilder2D()
+        builder = mesh_data.BoxMeshBuilder2D(tpe)
         surface_area = 1.0
     elif name == "box3d":
-        builder = mesh_data.BoxMeshBuilder3D()
+        builder = mesh_data.BoxMeshBuilder3D(tpe)
         surface_area = 1.0
     else:
         raise ValueError("unknown geometry name: %s" % name)
@@ -976,11 +985,11 @@ def rhs(t, w):
 # {{{ models: variable coefficient advection oversampling
 
 @pytest.mark.parametrize("order", [2, 3, 4])
-def test_improvement_quadrature(actx_factory, order):
+@pytest.mark.parametrize("tpe", [False, True])
+def test_improvement_quadrature(actx_factory, order, tpe):
     """Test whether quadrature improves things and converges"""
     from grudge.models.advection import VariableCoefficientAdvectionOperator
     from pytools.convergence import EOCRecorder
-    from meshmode.discretization.poly_element import QuadratureSimplexGroupFactory
     from meshmode.mesh import BTAG_ALL
 
     actx = actx_factory()
@@ -1002,23 +1011,30 @@ def conv_test(descr, use_quad):
         else:
             qtag = None
 
+        group_cls = TensorProductElementGroup if tpe else None
+        qfac = 2 if tpe else 4
         ns = [20, 25]
+        discr_order = order
         for n in ns:
             mesh = mgen.generate_regular_rect_mesh(
                 a=(-0.5,)*dims,
                 b=(0.5,)*dims,
                 nelements_per_axis=(n,)*dims,
-                order=order)
+                order=order, group_cls=group_cls)
 
             if use_quad:
                 discr_tag_to_group_factory = {
-                    qtag: QuadratureSimplexGroupFactory(order=4*order)
+                    dof_desc.DISCR_TAG_BASE:
+                    InterpolatoryEdgeClusteredGroupFactory(order),
+                    dof_desc.DISCR_TAG_QUAD:
+                    QuadratureGroupFactory(order=qfac*order)
                 }
+                discr_order = None
             else:
                 discr_tag_to_group_factory = {}
 
             dcoll = DiscretizationCollection(
-                actx, mesh, order=order,
+                actx, mesh, order=discr_order,
                 discr_tag_to_group_factory=discr_tag_to_group_factory
             )
 
@@ -1050,9 +1066,11 @@ def zero_inflow(dtag, t=0):
     eoc, errs = conv_test("no quadrature", False)
     q_eoc, q_errs = conv_test("with quadrature", True)
 
-    assert q_eoc > eoc
     assert (q_errs < errs).all()
     assert q_eoc > order - 0.1
+    # Fails for all tensor-product element types
+    assert q_eoc > eoc
+
 
 # }}}