Hi there,
I have a finite element theory question (so likely for @tjhei or @bangerth, but anyone with FE experience is welcome to chime in):
By default we often use a Gauss quadrature with FE_degree+1 for our assemblers, because that can accurately integrate polynomials of order 2*n-1 = 2*FE_degree+1, and can therefore accurately represent terms like phi_u * phi_u, which is of polynomial degree 2 * FE_degree. However we often use the same quadrature degrees in our postprocessors, where we simply want to integrate over the solution vector, which is of degree FE_degree. If I am not mistake in order to integrate this accurately we would only need a Gauss quadrature of order FE_degree/2 + 1. Is this correct? Because then we can reduce the quadrature order in many postprocessors (I am particularly looking at the depth average, because it is so expensive, and @maxrudolph needs it frequently in his models, every reduction in quadrature order would make it significantly cheaper).
Any advice is appreciated.
Rene
The way it is often taught is that a Gauss formula with q degrees is exact for
polynomials of degree 2q-1, and so to integrate finite element matrices like
the mass matrix (which has polynomial degree 2p), it is enough to use q=p+1.
In your case, you even have only polynomials of degree p, so q=roundup(p/2)
should be enough.
But this argument is not quite correct. It is only true for linear mappings
(on triangles, for example) because then the integrand really just remains a
polynomial of the same degree as on the reference cell. On the other hand, if
one is on quadrilaterals, or uses higher order mappings, then the actual
integrand one is integrating is no longer polynomial (but in reality some sort
of rational function). Consequently, Gauss quadrature of sufficiently high
order will no longer integrate exactly. Instead, what happens if you use
q=p+1 quadrature points is that the error one commits by replacing integration
by quadrature is of the same order as the finite element error, and so
sufficient to retain the overall convergence order. This is good enough. Using
a better quadrature formula would only reduce the constant in the error
(though likely not by very much in practice).
So your question is not so easy to answer. In practice, I suspect that if you
are integrating a quantity that is linear in the solution – say, the average
velocity/temperature/heat flux/… – then you are correct that a smaller
number of quadrature points will suffice. You won’t get the exact same answer,
but you shouldn’t get something that has a fundamentally different convergence
rate.
Best
W.
Thanks for the explanation!
I checked with the reduced quadrature order q=roundup((p+1)/2), and the difference in output is below 1e-3 for nearly all tests. There were some problems with exotic geometries (see https://github.com/geodynamics/aspect/pull/2910), but I suspect that is because of the geometries or model setup, not the postprocessing.
The depth averaging for SphericalShell and Box is nearly as accurate and significantly faster with the new version.
I checked with the reduced quadrature order |q=roundup((p+1)/2)|, and
the difference in output is below 1e-3 for nearly all tests.
Is this the relative difference? So is the difference in the third
significant digit of results? That really seems to be quite small, in
particular because almost all tests run on pretty coarse grids and the
error due to the discretization should still be quite large on these grids.
Yes that is relative difference on a coarse grid. Good enough I would say 