the code to generate non-negative weights for THIRD-order Gauss rule seems to be broken

[xxjaehoxx]

Dear all,

I tried setting the option allow_rules_with_negative_weights = FALSE for THIRD-order Gauss-like rule for my mesh with TET10 elements, and it seems to crash LAPACK. The actual error message that I get is: "the algorithm for computing the SVD failed to converge". Do you know what may be the issue here? Thank you!

Best, Mike

[jwpeterson]

In the case you mentioned, we fall back on a 2x2x2 conical product rule to avoid having negative weights. I wonder if there is some issue with the code that swaps between rules... if you can provide a simple test code that reproduces the issue, that would speed up my testing, otherwise I'll see if I can reproduce it.

[jwpeterson]

After reading your question a bit more carefully, it seems likely that you are producing a singular matrix by under-integrating. For example, the mass matrix for a Tet10 element using THIRD-order Gauss quadrature (negative weights or not) will be singular. This gist outputs:

n_qps = 8
mass_matrix = 0.00158519 0.000185185 6.66667e-05 -0.00017037 -0.00133333 -0.0022963 -0.00115556 -0.000681481 -0.0022963 -0.00223704 
0.000185185 0.00111111 0.000185185 0.000185185 -0.000740741 -0.000740741 -0.00259259 -0.00259259 -0.000740741 -0.00259259 
6.66667e-05 0.000185185 0.00134815 6.66667e-05 -0.0022963 -0.00133333 -0.000681481 -0.00271111 -0.0022963 -0.000681481 
-0.00017037 0.000185185 6.66667e-05 0.00158519 -0.0022963 -0.0022963 -0.00223704 -0.000681481 -0.00133333 -0.00115556 
-0.00133333 -0.000740741 -0.0022963 -0.0022963 0.0118519 0.00592593 0.00651852 0.00651852 0.00592593 0.00325926 
-0.0022963 -0.000740741 -0.00133333 -0.0022963 0.00592593 0.0118519 0.00651852 0.00325926 0.00592593 0.00651852 
-0.00115556 -0.00259259 -0.000681481 -0.00223704 0.00651852 0.00651852 0.0113778 0.00663704 0.00325926 0.00568889 
-0.000681481 -0.00259259 -0.00271111 -0.000681481 0.00651852 0.00325926 0.00663704 0.0104296 0.00651852 0.00663704 
-0.0022963 -0.000740741 -0.0022963 -0.00133333 0.00592593 0.00592593 0.00325926 0.00651852 0.0118519 0.00651852 
-0.00223704 -0.00259259 -0.000681481 -0.00115556 0.00325926 0.00651852 0.00568889 0.00663704 0.00651852 0.0113778 

Eigenvalue 0, (0.0416667, 0)
Eigenvalue 1, (0.00934363, 0)
Eigenvalue 2, (0.00881523, 0)
Eigenvalue 3, (0.00833333, 0)
Eigenvalue 4, (0.00146994, 0)
Eigenvalue 5, (0.00218477, 0)
Eigenvalue 6, (0.0025568, 0)
Eigenvalue 7, (5.90387e-19, 0)
Eigenvalue 8, (-1.96224e-18, 0)
Eigenvalue 9, (2.36944e-19, 0)

As you can see, there are three zero eigenvalues, so the mass matrix is not SPD.

Can you try using a higher-order quadrature and see if that fixes the issue that you observed? In the case of TET10, FIFTH order quadrature is sufficient to guarantee a non-singular mass matrix.

[boyceg]

@jwpeterson, this is not for constructing a mass matrix, just integrating a right-hand side. We are noticing some stability issues with THIRD order Gauss rules that do not seem to occur for SECOND or FIFTH order rules, and we are wondering if it is because the default THIRD order Gauss rule has negative weights. I believe that the only change that is made to generate this error is by changing allow_rules_with_negative_weights to false when setting up the quadrature rule. I think the code runs OK (albeit with compromised stability) if @xxjaehoxx leaves allow_rules_with_negative_weights at its default of true.

@xxjaehoxx, can you run in a debugger and get a stack trace for where this error is being emitted (presumably from the call to LAPACKgelss_ in DenseMatrix<T>::_svd_solve_lapack()) and/or try to create a simple standalone libMesh tester that can demonstrate the problem (e.g. I think it should happen just by trying to initialize the quadrature rule).

[jwpeterson]

We are noticing some stability issues

As I understand it, the main issue for rules with negative weights, especially if integrating "rough" data, is that even if the data is strictly positive, the quadrature rule can produce a non-positive integral. This could of course lead to numerical stability issues, depending on the application... I have never really seen this be an issue in practice, however, so it would definitely be interesting to see such a case.

Of course, if there is some library bug with swapping out quadrature rules on the fly as we do, that would be good to fix as well.

[jwpeterson]

The other thing that might be simple to try is changing your code to use a QConical rule at THIRD-order instead of QGauss. This would rule out the possibility of a library bug related to the allow_rules_with_negative_weights option.

[jwpeterson]

I also checked and we don't actually call svd_solve() anywhere in the library, so I guess the call that's failing must be in your own code somewhere? If you can describe generally how the matrix being SV decomposed is assembled, that might help me come up with some other debugging ideas...

[xxjaehoxx]

@jwpeterson sorry for the late reply. It took me some time to investigate what's going on. You're right. It seems like the call is failing on our side whenever we set allow_rules_with_negative_weights = FALSE. Thank you!