test_gradient: test warped mesh, overintegration

.github/workflows/ci.yml

@@ -12,9 +12,9 @@ jobs:
         name: Flake8
         runs-on: ubuntu-latest
         steps:
-        -   uses: actions/checkout@v3
+        -   uses: actions/checkout@v4
         -
-            uses: actions/setup-python@v4
+            uses: actions/setup-python@v5
             with:
                 # matches compat target in setup.py
                 python-version: '3.8'
@@ -27,7 +27,7 @@ jobs:
         name: Pylint
         runs-on: ubuntu-latest
         steps:
-        -   uses: actions/checkout@v3
+        -   uses: actions/checkout@v4
         -   name: "Main Script"
             run: |
                 echo "- matplotlib" >> .test-conda-env-py3.yml
@@ -46,7 +46,7 @@ jobs:
         name: Pytest on Py3
         runs-on: ubuntu-latest
         steps:
-        -   uses: actions/checkout@v3
+        -   uses: actions/checkout@v4
         -   name: "Main Script"
             run: |
                 curl -L -O https://gitlab.tiker.net/inducer/ci-support/raw/main/build-and-test-py-project-within-miniconda.sh
@@ -56,7 +56,7 @@ jobs:
         name: Examples on Py3
         runs-on: ubuntu-latest
         steps:
-        -   uses: actions/checkout@v3
+        -   uses: actions/checkout@v4
         -   name: "Main Script"
             run: |
                 curl -L -O https://tiker.net/ci-support-v0
@@ -78,7 +78,7 @@ jobs:
         name: Documentation
         runs-on: ubuntu-latest
         steps:
-        -   uses: actions/checkout@v3
+        -   uses: actions/checkout@v4
         -   name: "Main Script"
             run: |
                 echo "- matplotlib" >> .test-conda-env-py3.yml
@@ -100,7 +100,7 @@ jobs:
         name: Tests for downstream project ${{ matrix.downstream_project }}
         runs-on: ubuntu-latest
         steps:
-        -   uses: actions/checkout@v3
+        -   uses: actions/checkout@v4
         -   name: "Main Script"
             env:
                 DOWNSTREAM_PROJECT: ${{ matrix.downstream_project }}

grudge/discretization.py

@@ -35,6 +35,7 @@
 """
 
 from typing import Mapping, Optional, Union, TYPE_CHECKING, Any
+from meshmode.discretization.poly_element import ModalGroupFactory
 
 from pytools import memoize_method, single_valued
 
@@ -107,11 +108,8 @@ def _normalize_discr_tag_to_group_factory(
     assert discr_tag_to_group_factory is not None
 
     # Modal discr should always come from the base discretization
-    if DISCR_TAG_MODAL not in discr_tag_to_group_factory:
-        discr_tag_to_group_factory[DISCR_TAG_MODAL] = \
-            _generate_modal_group_factory(
-                discr_tag_to_group_factory[DISCR_TAG_BASE]
-            )
+    if DISCR_TAG_MODAL not in discr_tag_to_group_factory and order is not None:
+        discr_tag_to_group_factory[DISCR_TAG_MODAL] = ModalGroupFactory(order)
 
     return discr_tag_to_group_factory
 
@@ -391,15 +389,15 @@ def _base_to_geoderiv_connection(self, dd: DOFDesc) -> DiscretizationConnection:
         base_group_factory = self.group_factory_for_discretization_tag(
                 dd.discretization_tag)
 
-        def geo_group_factory(megrp, index):
+        def geo_group_factory(megrp):
             from modepy.shapes import Simplex
             from meshmode.discretization.poly_element import \
                     PolynomialEquidistantSimplexElementGroup
             if megrp.is_affine and issubclass(megrp._modepy_shape_cls, Simplex):
                 return PolynomialEquidistantSimplexElementGroup(
-                        megrp, order=0, index=index)
+                        megrp, order=0)
             else:
-                return base_group_factory(megrp, index)
+                return base_group_factory(megrp)
 
         from meshmode.discretization import Discretization
         geo_deriv_discr = Discretization(
@@ -779,30 +777,6 @@ def normal(self, dd):
     # }}}
 
 
-# {{{ modal group factory
-
-def _generate_modal_group_factory(nodal_group_factory):
-    from meshmode.discretization.poly_element import (
-        ModalSimplexGroupFactory,
-        ModalTensorProductGroupFactory
-    )
-    from meshmode.mesh import SimplexElementGroup, TensorProductElementGroup
-
-    order = nodal_group_factory.order
-    mesh_group_cls = nodal_group_factory.mesh_group_class
-
-    if mesh_group_cls is SimplexElementGroup:
-        return ModalSimplexGroupFactory(order=order)
-    elif mesh_group_cls is TensorProductElementGroup:
-        return ModalTensorProductGroupFactory(order=order)
-    else:
-        raise ValueError(
-            f"Unknown mesh element group: {mesh_group_cls}"
-        )
-
-# }}}
-
-
 # {{{ make_discretization_collection
 
 MeshOrDiscr = Union[Mesh, Discretization]

grudge/trace_pair.py

@@ -318,9 +318,10 @@ def interior_trace_pair(dcoll: DiscretizationCollection, vec) -> TracePair:
     return local_interior_trace_pair(dcoll, vec)
 
 
-def interior_trace_pairs(dcoll: DiscretizationCollection, vec, *,
-        comm_tag: Optional[Hashable] = None, tag: Hashable = None,
-        volume_dd: Optional[DOFDesc] = None) -> List[TracePair]:
+def interior_trace_pairs(
+        dcoll: DiscretizationCollection, vec, *,
+        comm_tag: Optional[Hashable] = None, volume_dd: Optional[DOFDesc] = None
+        ) -> List[TracePair]:
     r"""Return a :class:`list` of :class:`TracePair` objects
     defined on the interior faces of *dcoll* and any faces connected to a
     parallel boundary.
@@ -338,15 +339,6 @@ def interior_trace_pairs(dcoll: DiscretizationCollection, vec, *,
     :returns: a :class:`list` of :class:`TracePair` objects.
     """
 
-    if tag is not None:
-        warn("Specifying 'tag' is deprecated and will stop working in July of 2022. "
-                "Specify 'comm_tag' instead.", DeprecationWarning, stacklevel=2)
-        if comm_tag is not None:
-            raise TypeError("may only specify one of 'tag' and 'comm_tag'")
-        else:
-            comm_tag = tag
-    del tag
-
     if volume_dd is None:
         volume_dd = DD_VOLUME_ALL
 

test/test_modal_connections.py

@@ -23,12 +23,15 @@
 
 from grudge.array_context import PytestPyOpenCLArrayContextFactory
 from arraycontext import pytest_generate_tests_for_array_contexts
+
+from grudge.discretization import make_discretization_collection
 pytest_generate_tests = pytest_generate_tests_for_array_contexts(
         [PytestPyOpenCLArrayContextFactory])
 
 from meshmode.discretization.poly_element import (
     # Simplex group factories
     InterpolatoryQuadratureSimplexGroupFactory,
+    ModalGroupFactory,
     PolynomialWarpAndBlend2DRestrictingGroupFactory,
     PolynomialEquidistantSimplexGroupFactory,
     # Tensor product group factories
@@ -39,7 +42,6 @@
 from meshmode.dof_array import flat_norm
 import meshmode.mesh.generation as mgen
 
-from grudge import DiscretizationCollection
 import grudge.dof_desc as dof_desc
 
 import pytest
@@ -65,15 +67,16 @@ def f(x):
         group_cls=nodal_group_factory.mesh_group_class
     )
 
-    dcoll = DiscretizationCollection(
+    dcoll = make_discretization_collection(
         actx, mesh,
         discr_tag_to_group_factory={
-            dof_desc.DISCR_TAG_BASE: nodal_group_factory(order)
+            dof_desc.DISCR_TAG_BASE: nodal_group_factory(order),
+            dof_desc.DISCR_TAG_MODAL: ModalGroupFactory(order),
         }
     )
 
-    dd_modal = dof_desc.DD_VOLUME_MODAL
-    dd_volume = dof_desc.DD_VOLUME
+    dd_modal = dof_desc.DD_VOLUME_ALL_MODAL
+    dd_volume = dof_desc.DD_VOLUME_ALL
 
     x_nodal = actx.thaw(dcoll.discr_from_dd(dd_volume).nodes()[0])
     nodal_f = f(x_nodal)
@@ -105,17 +108,18 @@ def f(x):
         group_cls=QuadratureSimplexGroupFactory.mesh_group_class
     )
 
-    dcoll = DiscretizationCollection(
+    dcoll = make_discretization_collection(
         actx, mesh,
         discr_tag_to_group_factory={
             dof_desc.DISCR_TAG_BASE:
             PolynomialWarpAndBlend2DRestrictingGroupFactory(order),
-            dof_desc.DISCR_TAG_QUAD: QuadratureSimplexGroupFactory(2*order)
+            dof_desc.DISCR_TAG_QUAD: QuadratureSimplexGroupFactory(2*order),
+            dof_desc.DISCR_TAG_MODAL: ModalGroupFactory(order),
         }
     )
 
     # Use dof descriptors on the quadrature grid
-    dd_modal = dof_desc.DD_VOLUME_MODAL
+    dd_modal = dof_desc.DD_VOLUME_ALL_MODAL
     dd_quad = dof_desc.DOFDesc(dof_desc.DTAG_VOLUME_ALL,
                                dof_desc.DISCR_TAG_QUAD)
 

test/test_op.py

@@ -21,19 +21,26 @@
 """
 
 
+from meshmode.discretization.poly_element import (
+    InterpolatoryEdgeClusteredGroupFactory, QuadratureGroupFactory)
 import numpy as np
 
+from meshmode.mesh import BTAG_ALL
 import meshmode.mesh.generation as mgen
 
 from pytools.obj_array import make_obj_array
 
 from grudge import op, geometry as geo, DiscretizationCollection
-from grudge.dof_desc import DOFDesc
+from grudge.discretization import make_discretization_collection
+from grudge.dof_desc import (DISCR_TAG_BASE, DISCR_TAG_QUAD, DTAG_VOLUME_ALL,
+    DOFDesc, as_dofdesc)
 
 import pytest
 
 from grudge.array_context import PytestPyOpenCLArrayContextFactory
 from arraycontext import pytest_generate_tests_for_array_contexts
+
+from grudge.trace_pair import bv_trace_pair
 pytest_generate_tests = pytest_generate_tests_for_array_contexts(
         [PytestPyOpenCLArrayContextFactory])
 
@@ -44,37 +51,49 @@
 
 # {{{ gradient
 
-@pytest.mark.parametrize("form", ["strong", "weak"])
+@pytest.mark.parametrize("form", ["strong", "weak", "weak-overint"])
 @pytest.mark.parametrize("dim", [1, 2, 3])
 @pytest.mark.parametrize("order", [2, 3])
+@pytest.mark.parametrize("warp_mesh", [False, True])
 @pytest.mark.parametrize(("vectorize", "nested"), [
     (False, False),
     (True, False),
     (True, True)
     ])
 def test_gradient(actx_factory, form, dim, order, vectorize, nested,
-        visualize=False):
+                  warp_mesh, visualize=False):
     actx = actx_factory()
 
     from pytools.convergence import EOCRecorder
     eoc_rec = EOCRecorder()
 
-    for n in [4, 6, 8]:
-        mesh = mgen.generate_regular_rect_mesh(
-                a=(-1,)*dim, b=(1,)*dim,
-                nelements_per_axis=(n,)*dim)
+    for n in [8, 12, 16] if warp_mesh else [4, 6, 8]:
+        if warp_mesh:
+            if dim == 1:
+                pytest.skip("warped mesh in 1D not implemented")
+            mesh = mgen.generate_warped_rect_mesh(
+                          dim=dim, order=order, nelements_side=n)
+        else:
+            mesh = mgen.generate_regular_rect_mesh(
+                    a=(-1,)*dim, b=(1,)*dim,
+                    nelements_per_axis=(n,)*dim)
 
-        dcoll = DiscretizationCollection(actx, mesh, order=order)
+        dcoll = make_discretization_collection(
+                   actx, mesh,
+                   discr_tag_to_group_factory={
+                       DISCR_TAG_BASE: InterpolatoryEdgeClusteredGroupFactory(order),
+                       DISCR_TAG_QUAD: QuadratureGroupFactory(3 * order)
+                   })
 
         def f(x):
-            result = dcoll.zeros(actx) + 1
+            result = 1
             for i in range(dim-1):
                 result = result * actx.np.sin(np.pi*x[i])
             result = result * actx.np.cos(np.pi/2*x[dim-1])
             return result
 
         def grad_f(x):
-            result = make_obj_array([dcoll.zeros(actx) + 1 for _ in range(dim)])
+            result = make_obj_array([1 for _ in range(dim)])
             for i in range(dim-1):
                 for j in range(i):
                     result[i] = result[i] * actx.np.sin(np.pi*x[j])
@@ -87,16 +106,15 @@ def grad_f(x):
             result[dim-1] = result[dim-1] * (-np.pi/2*actx.np.sin(np.pi/2*x[dim-1]))
             return result
 
-        x = actx.thaw(dcoll.nodes())
-
-        if vectorize:
-            u = make_obj_array([(i+1)*f(x) for i in range(dim)])
-        else:
-            u = f(x)
+        def vectorize_if_requested(vec):
+            if vectorize:
+                return make_obj_array([(i+1)*vec for i in range(dim)])
+            else:
+                return vec
 
         def get_flux(u_tpair):
             dd = u_tpair.dd
-            dd_allfaces = dd.with_dtag("all_faces")
+            dd_allfaces = dd.with_domain_tag("all_faces")
             normal = geo.normal(actx, dcoll, dd)
             u_avg = u_tpair.avg
             if vectorize:
@@ -108,22 +126,45 @@ def get_flux(u_tpair):
                 flux = u_avg * normal
             return op.project(dcoll, dd, dd_allfaces, flux)
 
-        dd_allfaces = DOFDesc("all_faces")
+        x = actx.thaw(dcoll.nodes())
+        u = vectorize_if_requested(f(x))
+
+        bdry_dd_base = as_dofdesc(BTAG_ALL)
+        bdry_x = actx.thaw(dcoll.nodes(bdry_dd_base))
+        bdry_u = vectorize_if_requested(f(bdry_x))
 
         if form == "strong":
             grad_u = (
                 op.local_grad(dcoll, u, nested=nested)
                 # No flux terms because u doesn't have inter-el jumps
                 )
-        elif form == "weak":
+        elif form.startswith("weak"):
+            assert form in ["weak", "weak-overint"]
+            if "overint" in form:
+                quad_discr_tag = DISCR_TAG_QUAD
+            else:
+                quad_discr_tag = DISCR_TAG_BASE
+
+            allfaces_dd_base = as_dofdesc("all_faces", quad_discr_tag)
+            vol_dd_base = as_dofdesc(DTAG_VOLUME_ALL)
+            vol_dd_quad = vol_dd_base.with_discr_tag(quad_discr_tag)
+            bdry_dd_quad = bdry_dd_base.with_discr_tag(quad_discr_tag)
+            allfaces_dd_quad = allfaces_dd_base.with_discr_tag(quad_discr_tag)
+
             grad_u = op.inverse_mass(dcoll,
-                -op.weak_local_grad(dcoll, u, nested=nested)  # pylint: disable=E1130
+                -op.weak_local_grad(dcoll, vol_dd_quad,
+                        op.project(dcoll, vol_dd_base, vol_dd_quad, u),
+                        nested=nested)
                 +  # noqa: W504
                 op.face_mass(dcoll,
-                    dd_allfaces,
-                    # Note: no boundary flux terms here because u_ext == u_int == 0
-                    sum(get_flux(utpair)
-                        for utpair in op.interior_trace_pairs(dcoll, u))
+                    allfaces_dd_quad,
+                    sum(get_flux(
+                        op.project_tracepair(dcoll, allfaces_dd_quad, utpair))
+                        for utpair in op.interior_trace_pairs(
+                                      dcoll, u, volume_dd=vol_dd_base))
+                    + get_flux(
+                        op.project_tracepair(dcoll, bdry_dd_quad,
+                                   bv_trace_pair(dcoll, bdry_dd_base, u, bdry_u)))
                 )
             )
         else:
@@ -138,6 +179,9 @@ def get_flux(u_tpair):
             expected_grad_u = grad_f(x)
 
         if visualize:
+            # the code below does not handle the vectorized case
+            assert not vectorize
+
             from grudge.shortcuts import make_visualizer
             vis = make_visualizer(dcoll, vis_order=order if dim == 3 else dim+3)