From dff85cf1b0005537fa19e0f5d4334d53236f66cf Mon Sep 17 00:00:00 2001
From: Andreas Kloeckner <inform@tiker.net>
Date: Fri, 29 Nov 2024 15:44:45 -0600
Subject: [PATCH] Drop support for scalar-only broadcasting in dt_utils

---
 grudge/dt_utils.py | 59 +++++++++++++++-------------------------------
 1 file changed, 19 insertions(+), 40 deletions(-)

diff --git a/grudge/dt_utils.py b/grudge/dt_utils.py
index 36aac52bf..3587e371f 100644
--- a/grudge/dt_utils.py
+++ b/grudge/dt_utils.py
@@ -296,50 +296,29 @@ 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()
-                        },
-                        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, strict=True)))
-    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",
+    if not actx.supports_nonscalar_broadcasting:
+        raise ValueError("Array context does not allow advanced indexing. "
+                 "This is no longer supported.")
+
+    # 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]
-                    ) / afgrp.nunit_dofs,
-                    tagged=(FirstAxisIsElementsTag(),))
+                        face_ae_i.shape[-1])),
+                tagged=(FirstAxisIsElementsTag(),))
 
-                for vgrp, afgrp, face_ae_i in zip(volm_discr.groups,
-                                                  face_discr.groups,
-                                                  face_areas, strict=True)
-            )
-        )
+            for vgrp, face_ae_i in zip(volm_discr.groups, face_areas, strict=True)))
 
     return actx.freeze(
             actx.tag(NameHint(f"dt_geometric_{dd.as_identifier()}"),