Evolve TPE integration (#1043)

* Fix up TPE, and add some tests.
Co-authored-by: Matt Smith <[email protected]>

examples/gas-in-box.py

@@ -47,7 +47,7 @@
 from mirgecom.utils import force_evaluation
 from mirgecom.io import make_init_message
 
-from mirgecom.integrators import rk4_step
+from mirgecom.integrators import rk4_step, euler_step
 from mirgecom.steppers import advance_state
 from mirgecom.boundary import (
     AdiabaticSlipBoundary,
@@ -102,8 +102,8 @@ def main(actx_class, use_esdg=False, use_tpe=False,
          mech_name="uiuc_7sp", transport_type=0,
          use_av=0, use_limiter=False, order=1,
          nscale=1, npassive_species=0, map_mesh=False,
-         rotation_angle=0, add_pulse=False,
-         mesh_filename=None):
+         rotation_angle=0, add_pulse=False, nsteps=20,
+         mesh_filename=None, euler_timestepping=False):
     """Drive the example."""
     if casename is None:
         casename = "gas-in-box"
@@ -132,11 +132,11 @@ def main(actx_class, use_esdg=False, use_tpe=False,
         from leap.rk import RK4MethodBuilder
         timestepper = RK4MethodBuilder("state")
     else:
-        timestepper = rk4_step
-    n_steps = 20
+        timestepper = euler_step if euler_timestepping else rk4_step
+
     current_cfl = 1.0
     current_dt = 1e-6
-    t_final = current_dt * n_steps
+    t_final = current_dt * nsteps
     current_t = 0
     constant_cfl = False
     temperature_tolerance = 1e-2
@@ -193,6 +193,7 @@ def main(actx_class, use_esdg=False, use_tpe=False,
 
         local_mesh, global_nelements = distribute_mesh(comm, generate_mesh)
         local_nelements = local_mesh.nelements
+
         if dim is None:
             dim = local_mesh.ambient_dim
 
@@ -221,8 +222,7 @@ def add_wonk(x: np.ndarray) -> np.ndarray:
             theta = rotation_angle/180.0 * np.pi
             local_mesh = rotate_mesh_around_axis(local_mesh, theta=theta)
 
-    dcoll = create_discretization_collection(actx, local_mesh, order=order,
-                                             tensor_product_elements=use_tpe)
+    dcoll = create_discretization_collection(actx, local_mesh, order=order)
     nodes = actx.thaw(dcoll.nodes())
     ones = dcoll.zeros(actx) + 1.
 
@@ -751,6 +751,10 @@ def my_rhs(t, stepper_state):
         help="use leap timestepper")
     parser.add_argument("--esdg", action="store_true",
         help="use entropy-stable dg for inviscid terms.")
+    parser.add_argument("--euler-timestepping", action="store_true",
+                        help="use euler timestepping")
+    parser.add_argument("--nsteps", type=int, default=20,
+                        help="number of timesteps to take")
     parser.add_argument("--numpy", action="store_true",
         help="use numpy-based eager actx.")
     parser.add_argument("-d", "--dimension", type=int, choices=[1, 2, 3],
@@ -818,7 +822,7 @@ def my_rhs(t, stepper_state):
 
     main(actx_class, use_esdg=args.esdg, dim=args.dimension,
          use_overintegration=args.overintegration or args.esdg,
-         use_leap=args.leap, use_tpe=args.tpe,
+         use_leap=args.leap, use_tpe=args.tpe, nsteps=args.nsteps,
          casename=args.casename, rst_filename=rst_filename,
          periodic_mesh=args.periodic, use_mixture=args.mixture,
          multiple_boundaries=args.boundaries,
@@ -828,6 +832,6 @@ def my_rhs(t, stepper_state):
          use_navierstokes=args.navierstokes, npassive_species=args.species,
          nscale=args.weak_scale, mech_name=args.mechanism_name,
          map_mesh=args.wonky, rotation_angle=args.rotate, add_pulse=args.pulse,
-         mesh_filename=args.meshfile)
+         mesh_filename=args.meshfile, euler_timestepping=args.euler_timestepping)
 
 # vim: foldmethod=marker

mirgecom/discretization.py

@@ -49,35 +49,26 @@ def create_discretization_collection(actx, volume_meshes, order, *,
             DeprecationWarning, stacklevel=2)
 
     if tensor_product_elements:
-        warn("Overintegration is not supported for tensor product elements.")
+        warn(
+            "tensor_product_elements argument is not needed and will vanish soon.",
+            DeprecationWarning, stacklevel=2)
 
     from grudge.dof_desc import DISCR_TAG_BASE, DISCR_TAG_QUAD, DISCR_TAG_MODAL
     from grudge.discretization import make_discretization_collection
     from meshmode.discretization.poly_element import (
-        QuadratureSimplexGroupFactory,
-        PolynomialRecursiveNodesGroupFactory,
-        LegendreGaussLobattoTensorProductGroupFactory as Lgl,
+        InterpolatoryEdgeClusteredGroupFactory,
+        QuadratureGroupFactory,
         ModalGroupFactory
     )
 
     if quadrature_order < 0:
-        quadrature_order = 2*order+1
+        quadrature_order = 2*order + 1
 
-    if tensor_product_elements:
-        return make_discretization_collection(
-            actx, volume_meshes,
-            discr_tag_to_group_factory={
-                DISCR_TAG_BASE: Lgl(order),
-                DISCR_TAG_MODAL: ModalGroupFactory(order)
-            }
-        )
-    else:
-        return make_discretization_collection(
-            actx, volume_meshes,
-            discr_tag_to_group_factory={
-                DISCR_TAG_BASE: PolynomialRecursiveNodesGroupFactory(order=order,
-                                                                     family="lgl"),
-                DISCR_TAG_QUAD: QuadratureSimplexGroupFactory(quadrature_order),
-                DISCR_TAG_MODAL: ModalGroupFactory(order)
-            }
-        )
+    return make_discretization_collection(
+        actx, volume_meshes,
+        discr_tag_to_group_factory={
+            DISCR_TAG_BASE: InterpolatoryEdgeClusteredGroupFactory(order),
+            DISCR_TAG_MODAL: ModalGroupFactory(order),
+            DISCR_TAG_QUAD: QuadratureGroupFactory(quadrature_order)
+        }
+    )

test/test_operators.py

@@ -75,17 +75,17 @@ def _coord_test_func(dim, order=1):
 def _trig_test_func(dim):
     """Make trig test function.
 
-    1d: cos(2pi x)
-    2d: sin(2pi x)cos(2pi y)
-    3d: sin(2pi x)sin(2pi y)cos(2pi z)
+    1d: cos(pi x/2)
+    2d: sin(pi x/2)cos(pi y/2)
+    3d: sin(pi x/2)sin(pi y/2)cos(pi z/2)
     """
     sym_coords = prim.make_sym_vector("x", dim)
 
     sym_cos = pmbl.var("cos")
     sym_sin = pmbl.var("sin")
-    sym_result = sym_cos(2*np.pi*sym_coords[dim-1])
+    sym_result = sym_cos(np.pi*sym_coords[dim-1]/2)
     for i in range(dim-1):
-        sym_result = sym_result * sym_sin(2*np.pi*sym_coords[i])
+        sym_result = sym_result * sym_sin(np.pi*sym_coords[i]/2)
 
     return sym_result
 
@@ -133,8 +133,22 @@ def central_flux_boundary(actx, dcoll, soln_func, dd_bdry):
     return op.project(dcoll, bnd_tpair.dd, dd_allfaces, flux_weak)
 
 
-@pytest.mark.parametrize("tpe", [False, True])
-@pytest.mark.parametrize("dim", [1, 2, 3])
+@pytest.mark.parametrize(("dim", "mesh_name", "rot_axis", "wonk"),
+                         [
+                             (1, "tet_box1", None, False),
+                             (2, "tet_box2", None, False),
+                             (3, "tet_box3", None, False),
+                             (2, "hex_box2", None, False),
+                             (3, "hex_box3", None, False),
+                             (2, "tet_box2_rot", np.array([0, 0, 1]), False),
+                             (3, "tet_box3_rot1", np.array([0, 0, 1]), False),
+                             (3, "tet_box3_rot2", np.array([0, 1, 1]), False),
+                             (3, "tet_box3_rot3", np.array([1, 1, 1]), False),
+                             (2, "hex_box2_rot", np.array([0, 0, 1]), False),
+                             (3, "hex_box3_rot1", np.array([0, 0, 1]), False),
+                             (3, "hex_box3_rot2", np.array([0, 1, 1]), False),
+                             (3, "hex_box3_rot3", np.array([1, 1, 1]), False),
+                             ])
 @pytest.mark.parametrize("order", [1, 2, 3])
 @pytest.mark.parametrize("sym_test_func_factory", [
     partial(_coord_test_func, order=0),
@@ -144,7 +158,8 @@ def central_flux_boundary(actx, dcoll, soln_func, dd_bdry):
     _trig_test_func,
     _cv_test_func
 ])
-def test_grad_operator(actx_factory, tpe, dim, order, sym_test_func_factory):
+def test_grad_operator(actx_factory, dim, mesh_name, rot_axis, wonk,
+                       order, sym_test_func_factory):
     """Test the gradient operator for sanity.
 
     Check whether we get the right answers for gradients of analytic functions with
@@ -157,36 +172,67 @@ def test_grad_operator(actx_factory, tpe, dim, order, sym_test_func_factory):
     """
     import pyopencl as cl
     from grudge.array_context import PyOpenCLArrayContext
+    from meshmode.mesh.processing import rotate_mesh_around_axis
+    from grudge.dt_utils import h_max_from_volume
+    from mirgecom.simutil import componentwise_norms
+    from arraycontext import flatten
+    from mirgecom.operators import grad_operator
+    from meshmode.mesh.processing import map_mesh
+
+    def add_wonk(x: np.ndarray) -> np.ndarray:
+        wonk_field = np.empty_like(x)
+        if len(x) >= 2:
+            wonk_field[0] = (
+                1.5*x[0] + np.cos(x[0])
+                + 0.1*np.sin(10*x[1]))
+            wonk_field[1] = (
+                0.05*np.cos(10*x[0])
+                + 1.3*x[1] + np.sin(x[1]))
+        else:
+            wonk_field[0] = 1.5*x[0] + np.cos(x[0])
+
+        if len(x) >= 3:
+            wonk_field[2] = x[2] + np.sin(x[0] / 2) / 2
+
+        return wonk_field
+
+    tpe = mesh_name.startswith("hex_")
+    rotation_angle = 32.0
+    theta = rotation_angle/180.0 * np.pi
 
     # This comes from array_context
     actx = actx_factory()
 
     if tpe:  # TPE requires *grudge* array context, not array_context
-        if dim == 1:  # TPE does not support dim=1
-            pytest.skip()
         ctx = cl.create_some_context()
         queue = cl.CommandQueue(ctx)
         actx = PyOpenCLArrayContext(queue)
 
     sym_test_func = sym_test_func_factory(dim)
 
-    tol = 1e-10 if dim < 3 else 1e-9
+    tol = 1e-10 if dim < 3 else 2e-9
+
     from pytools.convergence import EOCRecorder
     eoc = EOCRecorder()
 
-    for nfac in [1, 2, 4]:
+    for nfac in [2, 4, 8]:
+
+        # Make non-uniform spacings
+        b = tuple((2 ** n) for n in range(dim))
 
-        mesh = get_box_mesh(dim, a=0, b=1, n=nfac*3, tensor_product_elements=tpe)
+        mesh = get_box_mesh(dim, a=0, b=b, n=nfac*3, tensor_product_elements=tpe)
+        if wonk:
+            mesh = map_mesh(mesh, add_wonk)
+        if rot_axis is not None:
+            mesh = rotate_mesh_around_axis(mesh, theta=theta, axis=rot_axis)
 
         logger.info(
             f"Number of {dim}d elements: {mesh.nelements}"
         )
 
-        dcoll = create_discretization_collection(actx, mesh, order=order,
-                                                 tensor_product_elements=tpe)
+        dcoll = create_discretization_collection(actx, mesh, order=order)
 
         # compute max element size
-        from grudge.dt_utils import h_max_from_volume
         h_max = h_max_from_volume(dcoll)
 
         def sym_eval(expr, x_vec):
@@ -208,14 +254,13 @@ def sym_eval(expr, x_vec):
         test_data = test_func(nodes)
         exact_grad = grad_test_func(nodes)
 
-        from mirgecom.simutil import componentwise_norms
-
-        from arraycontext import flatten
         err_scale = max(flatten(componentwise_norms(dcoll, exact_grad, np.inf),
                                 actx))
 
-        if err_scale <= 1e-16:
+        print(f"{err_scale=}")
+        if err_scale <= 1e-10:
             err_scale = 1
+            print(f"Rescaling: {err_scale=}")
 
         print(f"{test_data=}")
         print(f"{exact_grad=}")
@@ -225,7 +270,6 @@ def sym_eval(expr, x_vec):
         test_data_flux_bnd = _elbnd_flux(dcoll, int_flux, bnd_flux,
                                          test_data_int_tpair, boundaries)
 
-        from mirgecom.operators import grad_operator
         dd_vol = as_dofdesc("vol")
         dd_allfaces = as_dofdesc("all_faces")
         test_grad = grad_operator(dcoll, dd_vol, dd_allfaces,
@@ -236,6 +280,7 @@ def sym_eval(expr, x_vec):
             max(flatten(
                 componentwise_norms(dcoll, test_grad - exact_grad, np.inf),
                 actx) / err_scale)
+        print(f"{actx.to_numpy(h_max)=}\n{actx.to_numpy(grad_err)=}")
 
         eoc.add_data_point(actx.to_numpy(h_max), actx.to_numpy(grad_err))