diff --git a/docs/source/quadrature.rst b/docs/source/quadrature.rst index 02f0bb7ec6..c65dcddf2a 100644 --- a/docs/source/quadrature.rst +++ b/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,11 +30,28 @@ 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 @@ -42,9 +59,18 @@ For integrals that do not specify a quadrature degree, the default may be keyed 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 diff --git a/tests/firedrake/regression/test_quadrature.py b/tests/firedrake/regression/test_quadrature.py index 42e7b4423d..225a4b244d 100644 --- a/tests/firedrake/regression/test_quadrature.py +++ b/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"]) diff --git a/tsfc/kernel_interface/common.py b/tsfc/kernel_interface/common.py index aed238c64f..40d585786f 100644 --- a/tsfc/kernel_interface/common.py +++ b/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