Allow passing a maximum degree quadrature to a measure

docs/source/quadrature.rst

@@ -15,9 +15,9 @@ this estimate might be quite large, and a warning like this one will be raised:
 
    tsfc:WARNING Estimated quadrature degree 13 more than tenfold greater than any argument/coefficient degree (max 1)
 
-For integrals with very complicated nonlinearities, the estimated quadrature
-degree might be in the hundreds or thousands, rendering the integration
-prohibitively expensive, or leading to segfaults.
+For integrals with very complicated or non-polynomial nonlinearities, the
+estimated quadrature degree might be in the hundreds or thousands, rendering
+the integration prohibitively expensive, or leading to segfaults.
 
 Specifying the quadrature rule in the variational formulation
 -------------------------------------------------------------
@@ -30,21 +30,47 @@ quadrature degree can be prescribed on each integral :py:class:`~ufl.measure.Mea
    inner(sin(u)**4, v) * dx(degree=4)
 
 Setting ``degree=4`` means that the quadrature rule will be exact only for integrands
-of total polynomial degree up to 4. This, of course, will introduce a greater numerical error than the default.
+of total polynomial degree up to 4. This, of course, will introduce a greater numerical
+error than the default.
 
-For integrals that do not specify a quadrature degree, the default may be keyed as
-``"quadrature_degree"`` in the ``form_compiler_parameters`` dictionary passed on to
-:py:func:`~.solve`, :py:func:`~.project`, or :py:class:`~.NonlinearVariationalProblem`.
+Rather than enforcing a fixed quadrature degree, it is also possible to set
+a maximum allowable degree. This value will only be used if UFL estimates a larger
+degree.
+The maximum allowable degree can be set for a particular integral by
+adding the ``"max_quadrature_degree"`` entry to the ``metadata`` of the ``Measure``:
+
+.. code-block:: python3
+
+   inner(sin(u)**4, v) * dx(metadata={"max_quadrature_degree": 4})
+
+The ``metadata`` can also be used to set a fixed quadrature degree by adding
+a ``"quadrature_degree"`` entry to the dictionary. However, because setting a
+fixed degree is quite common the ``degree`` keyword argument shown above is
+provided as a convenient shorthand.
+
+For integrals that do not specify a fixed or maximum quadrature degree, a default value
+may be keyed as the ``"quadrature_degree"`` or ``"max_quadrature_degree"`` entry
+respectively in the ``form_compiler_parameters`` dictionary passed on to :py:func:`~.solve`,
+:py:func:`~.project`, :py:class:`~.NonlinearVariationalProblem`, or :py:func:`~.assemble`.
 
 .. code-block:: python3
 
    F = inner(grad(u), grad(v))*dx(degree=0) + inner(exp(u), v)*dx - inner(1, v)*dx
 
    solve(F == 0, u, form_compiler_parameters={"quadrature_degree": 4})
 
-In the example above, only the integrals with unspecified quadrature degree
-will be computed on a quadrature rule that exactly integrates polynomials of
-the degree set in ``form_compiler_parameters``.
+   assemble(F, form_compiler_parameters={"max_quadrature_degree": 4})
+
+In the example above:
+
+* In the ``solve`` call, the integrals with unspecified quadrature degree
+  will be computed on a quadrature rule that exactly integrates polynomials of
+  the degree set in ``form_compiler_parameters``.
+
+* In the ``assemble`` call, the integrals with unspecified quadrature degree will
+  be computed with the ``"max_quadrature_degree"`` set in the ``"form_compiler_parameters"``
+  if the estimated quadrature degree is greater than 4, otherwise they will be computed
+  with the estimated quadrature degree.
 
 Another way to specify the quadrature rule is through the ``scheme`` keyword. This could be
 either a :py:class:`~finat.quadrature.QuadratureRule`, or a string. Supported string values

tests/firedrake/regression/test_quadrature.py

@@ -20,6 +20,25 @@ def test_hand_specified_quadrature(mesh):
     assert not np.allclose(a_q0.dat.data, a_q2.dat.data)
 
 
+def test_hand_specified_max_quadrature():
+    mesh = UnitIntervalMesh(1)
+
+    x, = SpatialCoordinate(mesh)
+    a = (x**4)*dx
+
+    # These should be the same because we only need degree=4 for exact integration.
+    x4 = assemble(a)
+    x4_maxquad5 = assemble(a, form_compiler_parameters={"max_quadrature_degree": 5})
+
+    assert np.isclose(x4, x4_maxquad5)
+
+    # These should be the same because degree=2 will limit the quadrature
+    x4_quad2 = assemble(a, form_compiler_parameters={"quadrature_degree": 2})
+    x4_maxquad2 = assemble(a, form_compiler_parameters={"max_quadrature_degree": 2})
+
+    assert np.isclose(x4_quad2, x4_maxquad2)
+
+
 @pytest.mark.parametrize("diagonal", [False, True])
 @pytest.mark.parametrize("mat_type", ["matfree", "aij"])
 @pytest.mark.parametrize("family", ["Quadrature", "Boundary Quadrature"])

tsfc/kernel_interface/common.py

@@ -320,11 +320,20 @@ def set_quad_rule(params, cell, integral_type, functions):
                         if e.family() in {"Quadrature", "Boundary Quadrature"})
         if len(quad_data) == 0:
             quadrature_degree = params["estimated_polynomial_degree"]
-            if all((asarray(quadrature_degree) > 10 * asarray(e.degree())).all() for e in elements):
-                logger.warning("Estimated quadrature degree %s more "
-                               "than tenfold greater than any "
-                               "argument/coefficient degree (max %s)",
-                               quadrature_degree, max_degree([e.degree() for e in elements]))
+            if "max_quadrature_degree" in params:
+                max_allowed_degree = params["max_quadrature_degree"]
+                if quadrature_degree > max_allowed_degree:
+                    logger.info("Estimated quadrature degree %s greater "
+                                "than maximum allowed degree %s. "
+                                "Using maximum degree %s instead.",
+                                quadrature_degree, max_allowed_degree, max_allowed_degree)
+                    quadrature_degree = max_allowed_degree
+            else:
+                if all((asarray(quadrature_degree) > 10 * asarray(e.degree())).all() for e in elements):
+                    logger.warning("Estimated quadrature degree %s more "
+                                   "than tenfold greater than any "
+                                   "argument/coefficient degree (max %s)",
+                                   quadrature_degree, max_degree([e.degree() for e in elements]))
         else:
             try:
                 (quadrature_degree, quad_rule), = quad_data