Rule of thumb for quadrature degree

variational formulation
1

Hello everyone,

I have successfully implemented LC-PML in FEniCSx, but I need some help fine-tuning a few details. The equation being solved essentially becomes an anisotropic version of the Helmholtz equation:

\int_{\Omega}\nabla w \cdot \nabla p d\Omega - \int_{\Omega}k^{2}wp d\Omega + \int_{\Omega_{PML}}\Lambda_{PML} \nabla w \cdot \nabla p d\Omega_{PML} - \int_{\Omega_{PML}}\alpha_{PML}k^{2}wp d\Omega_{PML} = \int_{\Omega}j\omega\rho_{0} w q_{m} d\Omega +\oint_{\partial\Omega}w \frac{\partial p}{\partial n} dS

Here, \Lambda_{PML} and \alpha_{PML} define the behavior of the PML.

Since these functions can be quite complex and differ significantly from polynomials, ffcx is estimating a very high quadrature degree (between 19 and 25).

I am aware that I can manually set the quadrature degree using the following:

dx = Measure("dx", domain=msh, metadata={"quadrature_degree": desired_degree})

However, I would appreciate some advice:

  1. Is it common for complicated weak forms to result in a high estimated quadrature degree?
  2. Can I use a rule of thumb to determine the appropriate quadrature degree based on the element order? Should it be more than twice the element order, considering the presence of the \Lambda_{PML} tensor?

While I plan to conduct a sensitivity analysis (which I will do), your advice would help put my mind at ease.
Thank you!

2

There are a few things to consider here:

  1. is your element type triangle/tetrahedra or quadrilateral/hexahedral
    1b) are the mesh cells affine?

If no on 1b the degree is usually alot higher than in the affine case.

Consider a general affine simplex with a test/trial function of order Q. Then a mass matrix has order 2Q, while a stiffness matrix has order 2*(Q+1).

  1. How does the PML functions look like? In general you can use
    How to obtain the Gauss quadrature degree estimated by the form compiler - #2 by dokken
    to check what the estimated degree of various forms are.

Without an example on a square or cube with the forms (ignoring the subdomains of pml, just using regular dx), it is hard to give further pointers.

3
  1. I am using linear/quadratic tetrahedra
    1b. what is an affine mesh? Sorry if it should be obvious.
    The model is based on a complex-curvilinear stretch of the real coordinates, which is taken into account in the tensor \Lambda_{PML} = det(J_{PML})J_{PML}^{-1} J_{PML}^{-T}, and the coefficient \alpha_{PML} = det(J_{PML}) , where J_{PML} is the Jacobian matrix of the complex-curvilinear transformation.
    I am planning to make soon a detailed description of what I did, but everithing is based on the following:

https://www.researchgate.net/publication/345315346_An_automatic_PML_for_acoustic_finite_element_simulations_in_convex_domains_of_general_shape

https://www.sciencedirect.com/science/article/abs/pii/S0010465523000863

  1. the output of estimate_total_polynomial_degree() gives me 19, that results in the warning:
  WARNING: The number of integration points for each cell will be: 1000
           Consider using the option 'quadrature_degree' to reduce the number of points

The thing is: if I don’t reduce the quadrature degree, I get a huge computation time. Do I need such a giant amount of integration points?

Thank you very much

4

An affine mesh is a mesh Which has a linear operation between the reference and physical element.

First order triangles and tetrahedra are affine.
In general quads and hexes are not (only if they are parallellograms or a parallelepiped).
If they are non-affine, the jacobian is a polynomial of some order (in 3D i think its about 4).

Higher order (curved cells) are also non-affine.

I would use the Xiao-Gimbutas quadrature rule:

Which can be accessed through Basix and used as a custom quadrature rule.

Ref

5

It seems that this is defined only for triangles up to m=19, while for tetra it goes up to m=15.

From the code it seems that for higher values it switches to gauss-jacobi. Is it correct?

Anyway, I show you a sensitivity analysis on the SPL spectrum in a single point, with different quadrature degrees (ranging from 2 to 19). Of course everything else (including the order of the tetra) is the same. The result is not varying. Is it right? The only thing that doesn’t remain the same is the computation time, that goes from 2 minutes to 20 minutes.

comparison_quadrature_degrees
comparison_quadrature_degrees927×642 36.9 KB

6

Sorry, yes the upper limit is 15, but at least it is alot less points than in the case of use GLL.

As you haven’t provided the ufl formulation of the PML I cannot tell you why ufl estimates the quadrature degree to be this high.

Looking at your graph, indeed, it seems like ufl is over-estimating the quadrature degree (compared to what your problem needs).

7

The ufl formulation is identical to the formulas posted above:

it is:

a = inner(grad(u), grad(v)) * dx(1) - k0**2 * inner(u, v) * dx(1)  \
        + inner(LAMBDA_PML*grad(u), grad(v)) * dx(2) - detJ * k0**2 * inner(u, v) * dx(2)
L = inner(f, v) * dx(1)

starting from quantities that have been computed in gmsh, that are:

k2, k3 = principal curvatures of the PML interface
t2, t3 = principal directions of the PML interface
n      = normal vector of the PML interface
phi_domain = normalized distance from the PML interface 
d_PML = total thickness of PML layer

I get to LAMBDA_PML and detJ:

        C = 15*k0

        N_stretching = 3
        sigma = C*phi_domain**N_stretching
        f_csi = d_pml*C/(N_stretching+1)*phi_domain**(N_stretching+1)

        csi = phi_domain*d_pml

        h2 = 1 + k2*csi
        h3 = 1 + k3*csi

        s1 = 1 + 1/(1j*k0)*sigma
        s2 = 1 + k2 * f_csi/(1j*k0*h2) 
        s3 = 1 + k3 * f_csi/(1j*k0*h3) 

        nnT = as_matrix(np.outer(n,n))
        t2t2T = as_matrix(np.outer(t2,t2))
        t3t3T = as_matrix(np.outer(t3,t3))

        detJ = s1*s2*s3
        LAMBDA_PML = s2*s3/s1 * nnT + s1*s3/s2 * t2t2T + s1*s2/s3 * t3t3T

I am not posting the entire MWE because it is huge and not easy to read by human beings. I am working on making it cristal clear! Sorry