proposed changes to production

grudge/array_context.py

@@ -303,19 +303,21 @@ def _dag_to_compiled_func(self, dict_of_named_arrays,
 #                    dict_of_named_arrays)
 #            else:
 #                raise
-        distributed_partition = pt.find_distributed_partition(
-            # pylint-ignore-reason:
-            # '_BasePytatoArrayContext' has no
-            # 'mpi_communicator' member
-            # pylint: disable=no-member
-            self.actx.mpi_communicator, dict_of_named_arrays)
-
-        if __debug__:
-            # pylint-ignore-reason:
-            # '_BasePytatoArrayContext' has no 'mpi_communicator' member
-            pt.verify_distributed_partition(
-                self.actx.mpi_communicator,  # pylint: disable=no-member
-                distributed_partition)
+
+        with ProcessLogger(logger, "pt.find_distributed_partition"):
+            distributed_partition = pt.find_distributed_partition(
+                # pylint-ignore-reason:
+                # '_BasePytatoArrayContext' has no
+                # 'mpi_communicator' member
+                # pylint: disable=no-member
+                self.actx.mpi_communicator, dict_of_named_arrays)
+
+            if __debug__:
+                # pylint-ignore-reason:
+                # '_BasePytatoArrayContext' has no 'mpi_communicator' member
+                pt.verify_distributed_partition(
+                    self.actx.mpi_communicator,  # pylint: disable=no-member
+                    distributed_partition)
 
         self.actx._compile_trace_callback(self.f, "post_find_distributed_partition",
                 distributed_partition)

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

test/mesh_data.py

@@ -117,43 +117,38 @@ class _BoxMeshBuilderBase(MeshBuilder):
     group_cls = None
     a = (-0.5, -0.5, -0.5)
     b = (+0.5, +0.5, +0.5)
+    tpe: bool = False
+
+    def __init__(self, tpe=False, a=(-0.5, -0.5, -0.5),
+                 b=(0.5, 0.5, 0.5)):
+        self.tpe = tpe
+        self.a = a
+        self.b = b
 
     def get_mesh(self, resolution, mesh_order=None):
         if mesh_order is None:
             mesh_order = self.mesh_order
         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,
-                group_cls=self.group_cls,
-                order=mesh_order)
+                group_cls=group_cls, order=mesh_order)
 
 
 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 = 3
 
-    def __init__(self, tpe=False):
-        if tpe:
-            self.group_cls = TensorProductElementGroup
-
 
 class WarpedRectMeshBuilder(MeshBuilder):
     resolutions: ClassVar[Sequence[Hashable]] = [4, 6, 8]

test/test_dt_utils.py

@@ -107,12 +107,12 @@ def test_geometric_factors_regular_refinement(actx_factory, name, tpe):
     assert np.all(np.isclose(ratios, 2))
 
     # Make sure it works with empty meshes
-    if not tpe:
-        mesh = builder.get_mesh(0, order)
-        dcoll = make_discretization_collection(
-            actx, mesh, order=order,
-            discr_tag_to_group_factory=dtag_to_grp_fac)
-        factors = actx.thaw(dt_geometric_factors(dcoll))  # noqa: F841
+    # if not tpe:
+    #    mesh = builder.get_mesh(0, order)
+    #    dcoll = make_discretization_collection(
+    #        actx, mesh, order=order,
+    #        discr_tag_to_group_factory=dtag_to_grp_fac)
+    #    factors = actx.thaw(dt_geometric_factors(dcoll))  # noqa: F841
 
 
 @pytest.mark.parametrize("name", ["interval", "box2d", "box3d"])

test/test_grudge.py

@@ -158,17 +158,21 @@ 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"
+    "box2d-tpe", "box3d-tpe", "box2d", "box3d", "2-1-ellipse", "spheroid",
     ])
-def test_mass_surface_area(actx_factory, tpe, name):
+@pytest.mark.parametrize("oi", [False, True])
+def test_mass_surface_area(actx_factory, name, oi):
+    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)
+    vol_dd_quad = vol_dd_base.with_discr_tag(qtag)
 
     # {{{ cases
 
     order = 4
-
+    tpe = name.endswith("-tpe")
     if name == "2-1-ellipse":
         if tpe:
             pytest.skip()
@@ -179,11 +183,11 @@ def test_mass_surface_area(actx_factory, tpe, name):
             pytest.skip()
         builder = mesh_data.SpheroidMeshBuilder()
         surface_area = _spheroid_surface_area(builder.radius, builder.aspect_ratio)
-    elif name == "box2d":
-        builder = mesh_data.BoxMeshBuilder2D(tpe)
+    elif name.startswith("box2d"):
+        builder = mesh_data.BoxMeshBuilder2D(tpe=tpe)
         surface_area = 1.0
-    elif name == "box3d":
-        builder = mesh_data.BoxMeshBuilder3D(tpe)
+    elif name.startswith("box3d"):
+        builder = mesh_data.BoxMeshBuilder3D(tpe=tpe)
         surface_area = 1.0
     else:
         raise ValueError(f"unknown geometry name: {name}")
@@ -198,17 +202,24 @@ def test_mass_surface_area(actx_factory, tpe, 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)
+        nodes = actx.thaw(volume_discr.nodes())
 
         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))
+        print(f"approx_surface_area({name})={approx_surface_area}")
 
         logger.info(
             f"surface: got {approx_surface_area:.5e} / expected {surface_area:.5e}")  # noqa: G004
@@ -229,6 +240,154 @@ def test_mass_surface_area(actx_factory, tpe, name):
 
     assert eoc.max_error() < 3e-13 or eoc.order_estimate() > order
 
+
+@pytest.mark.parametrize("name", [
+    "interval", "box2d", "box2d-tpe", "box3d", "box3d-tpe"
+    ])
+# @pytest.mark.parametrize("discr_order", [1, 2, 3, 4, 5])
+def test_correctness_of_quadrature(actx_factory, name):
+    from grudge.dof_desc import DISCR_TAG_BASE, DISCR_TAG_QUAD, as_dofdesc
+    actx = actx_factory()
+    vol_dd_base = as_dofdesc(dof_desc.DTAG_VOLUME_ALL)
+    vol_dd_quad = vol_dd_base.with_discr_tag(DISCR_TAG_QUAD)
+
+    # {{{ cases
+
+    tol = 5e-13
+    # discr_order = 1
+    dim = None
+    mesh_order = 1
+
+    tpe = name.endswith("-tpe")
+    if name.startswith("box2d"):
+        builder = mesh_data.BoxMeshBuilder2D(tpe=tpe,
+                                             a=(0, 0),
+                                             b=(2.0, 2.0))
+        dim = 2
+    elif name.startswith("box3d"):
+        builder = mesh_data.BoxMeshBuilder3D(tpe=tpe,
+                                             a=(0, 0, 0),
+                                             b=(2.0, 2.0, 2.0))
+        dim = 3
+    elif name == "interval":
+        builder = mesh_data.BoxMeshBuilder1D(tpe=False, a=(0.0,),
+                                             b=(2.0,))
+        dim = 1.0
+    else:
+        raise ValueError(f"unknown geometry name: {name}")
+    exact_volume = 2.0**dim
+    print(f"Domain: {name} ({dim}d), {exact_volume=}")
+    print(f"======================================================")
+
+    # }}}
+
+    # {{{ convergence
+
+    from pytools.convergence import EOCRecorder
+    for discr_order in range(1, 8):
+        print(f"   {discr_order=}")
+        print("   --------------------")
+        report_discr = True
+        for field_order in range(1, max(2*discr_order+1, 8)):
+            report_field_order = True
+            eoc_base = EOCRecorder()
+            eoc_quad = EOCRecorder()
+            for resolution in builder.resolutions:
+                mesh = builder.get_mesh(resolution, mesh_order)
+                dcoll = make_discretization_collection(
+                    actx, mesh,
+                    discr_tag_to_group_factory={
+                        DISCR_TAG_BASE:
+                        InterpolatoryEdgeClusteredGroupFactory(discr_order),
+                        DISCR_TAG_QUAD: QuadratureGroupFactory(discr_order)
+                    })
+                vol_discr_base = dcoll.discr_from_dd(vol_dd_base)
+                vol_discr_quad = dcoll.discr_from_dd(vol_dd_quad)
+                if report_discr:
+                    nelem = vol_discr_base.mesh.nelements
+                    ndofs_base = vol_discr_base.ndofs
+                    ndofs_quad = vol_discr_quad.ndofs
+                    dofs_per_el_base = ndofs_base/nelem
+                    dofs_per_el_quad = ndofs_quad/nelem
+                    print(f"   - {dofs_per_el_base=}, {dofs_per_el_quad=}")
+                    report_discr = False
+                if report_field_order:
+                    print(f"      - {field_order=}")
+                    print("      - - - - - - - - - -")
+                    report_field_order = False
+                nodes_base = actx.thaw(vol_discr_base.nodes())
+                nodes_quad = actx.thaw(vol_discr_quad.nodes())
+                ones_base = 0*nodes_base[0] + 1
+                ones_quad = 0*nodes_quad[0] + 1
+
+                approx_vol_base = \
+                    actx.to_numpy(op.integral(dcoll, vol_dd_base, ones_base))
+                approx_vol_quad = \
+                    actx.to_numpy(op.integral(dcoll, vol_dd_quad, ones_quad))
+                err_vol_base = abs(approx_vol_base - exact_volume)/exact_volume
+                err_vol_quad = abs(approx_vol_quad - exact_volume)/exact_volume
+
+                logger.info(
+                    f"Name: {name} ({dim}d)\n"
+                    f"Exact volume: {exact_volume}\n"
+                    f"volume base: got {approx_vol_base:.12e}\n"  # noqa: G004
+                    f"volume quad: got {approx_vol_quad:.12e}\n")  # noqa: G004
+
+                # Quadrature should get exact volume for all discr (p >= 1)
+                assert err_vol_base < tol
+                assert err_vol_quad < tol
+
+                field_base = nodes_base[0]**field_order
+                field_quad = nodes_quad[0]**field_order
+                if dim > 1:
+                    field_base = sum(nodes_base[i]**field_order
+                                     for i in range(dim))
+                    field_quad = sum(nodes_quad[i]**field_order
+                                     for i in range(dim))
+                ofac = 1.0/2**field_order
+                field_base = ofac * field_base * (field_order + 1.0)
+                field_quad = ofac * field_quad * (field_order + 1.0)
+                exact_integral = dim*exact_volume
+
+                integral_base =  \
+                    actx.to_numpy(op.integral(dcoll, vol_dd_base, field_base))
+                integral_quad =  \
+                    actx.to_numpy(op.integral(dcoll, vol_dd_quad, field_quad))
+                err_base = \
+                    abs(integral_base - exact_integral)/exact_integral
+                err_quad = \
+                    abs(integral_quad - exact_integral)/exact_integral
+
+                if field_order <= discr_order:
+                    assert err_base < tol
+                    assert err_quad < tol
+
+                # compute max element size
+                from grudge.dt_utils import h_max_from_volume
+                h_max = h_max_from_volume(dcoll)
+
+                eoc_base.add_data_point(actx.to_numpy(h_max), err_base)
+                eoc_quad.add_data_point(actx.to_numpy(h_max), err_quad)
+
+            logger.info("volume error(base)\n%s", str(eoc_base))
+            logger.info("volume error(quad)\n%s", str(eoc_quad))
+            print("---- base -----")
+            print(f"{eoc_base.pretty_print()}")
+            print("---- quad -----")
+            print(f"{eoc_quad.pretty_print()}")
+
+            # Sanity check here: *must* be exact if discr_order is sufficient
+            if discr_order >= field_order:
+                assert eoc_base.max_error() < tol
+                assert eoc_quad.max_error() < tol
+            else:
+                if eoc_base.max_error() > tol:  # *can* be exact(ish) otherwise
+                    assert eoc_base.order_estimate() > discr_order
+                if eoc_quad.max_error() > tol:
+                    assert eoc_quad.order_estimate() > discr_order
+
+        print("-=-=-=-=-=-=-=-=-=-=-=-=-=--=-=-=-")
+    print("==============================")
 # }}}