Pardiso with NURBFEAP

[resammc]

Dear Prof. Taylor/Feap Admin,

Is there any special care which must be taken when using Pardiso solver with nurbfeap (or even feap itself)?

The problem is, when I use FEAP's standard solver, the solution converges, but using Pardiso it converges strangely and for some cases it doesn't converge (NaN or Inf residual after one iteration...).

The problem is not with FEAP's standard elements, but with my user elements, so I thought maybe I'm missing something here.

OS: Ubuntu 16.04
Compiler+Math Library: Intel® Parallel Studio XE Cluster Edition for Linux* 2016 Update 3

the problem occurs with both TANG and UTANG.

Any help regarding this would be greatly appreciated.

[Prof. S. Govindjee]

Pardiso should not care about this.  There are two things to check however
(1) Have a look at the matrix itself when just using FEAP's standard solver and compare to what is being passed to Pardiso
(2) Make sure there are no conflicts with user macro files when adding in Pardiso with those in nurbFEAP

[resammc]

Dear Prof. Govindjee,

Thank you very much for your kind response and the useful suggestions. I investigated more on the generated tangent arrays.

when I look at the TANG1 array (via SHOW command), size of this array is different with different solvers (pardiso uses a compressed form) (unfortunately, it is not possible to use OUTP,TANG command to extract Pardiso Tangent matrix) and a one-by-one comparison is not possible, however:

STANDARD FEAP SOLVER, TANG1 ARRAY (SIZE: 2588)

FIRST 10 COMPONENTS:

Code: [Select]

    Matrix: TANG1                         
    row/col     1           2           3           4           5           6
       1  1.7603E-05  9.7348E-06  1.3176E-05  1.0315E-05  1.3295E-05  1.1151E-05
       2  1.3091E-05  1.2249E-05  1.2050E-05  1.3407E-05
MKL PARDISO SOLVER, TANG1 ARRAY (SIZE: 1099)

FIRST 10 COMPONENTS:

Code: [Select]

    Matrix: TANG1                         
    row/col     1           2           3           4           5           6
       1  5.6808E+04  1.8702E+04  9.7255E+01 -5.1553E+04 -2.3957E+04 -1.7072E+04
       2  1.0888E+05  1.0071E+03  1.1376E+04  6.2694E+03
as you see, the tangent generated by Pardiso looks like an unfactorized tangent.

regarding the conflicts with NURBFEAP, I'm using umacr20 which is not reserved in nurbfeap.

P.S: there is no problem with the residual matrix + even TANG,NUME (numerical tangent) produces two different arrays for different solvers.

[FEAP_Admin]

If you use TANG,,-1 it will produce only the unfactored tangent.  Please try that.
And yes the storage format is different when you used padisio versus FEAP's profile solver.
But FEAP also has built-in solvers that use the pardisio format (I think).  Try using the
command DIREct,SPARse before generating the tangent.

[Prof. R.L. Taylor]

I agree with Professor Govindjee that there should be no differences.  A couple of queries,

Is your matrix symmetric or unsymmetric?

Are you sure the element produces the same element matrix for both forms?  I am sure there a re some flags that are set  for nurbs elements -- check the elements we produce with yours to see that they are all defined.

Follow the assembly process to see that pardiso is assmbling with the sparse format.

If your element is not too complicated maybe oosting or email to feap account will help us debug what is happening.

We will work on producing an output for sparse form so we can check the produced matrices are the same.

[Prof. S. Govindjee]

One other debugging mechanism could be

TANG,,-1
OUTP,TANG

or

TANG,,-1
OUTP,UTAN

depending on your problem.  Looking at program/pouta.f this should work for the sparse as well as profil storage methods.

[Prof. R.L. Taylor]

The attached file will add an option for Pardiso solver to ouput UNFACTORED matrix (i.e., you must user: TANG,,-1).  To output the file use the solution command: PARDiso OUTPut

[resammc]

Dear Prof. Taylor / Prof. Govindjee / FEAP Admin

Thank you very much for all the suggestions. I will go through them all and will report back here.

[resammc]

I went through some tests.

First, I used the output commands to print out the non-symmetric tangent (in my problem, the matrix is non-symmetric).

Investigated cases:

FEAP STANDARD - FACTORED: (UTANG,,1 -> OUTP UTAN)
FEAP STANDARD - UNFACTORED: (UTANG,,-1 -> OUTP UTAN)
PARDISO - FACTORED: (PARD -> UTANG,,1  -> PARD OUTP)
PARDISO - UNFACTORED: (PARD -> UTANG,,-1 -> PARD OUTP)

The outputs were not in a same order, so I tried to sort them. Here are the results (only some of the components are shown here, kindly look at the attachments for complete data):

For diagonal components:

                 FEAP_FACT   FEAP_UNFACT   PARD_FACT   PARD_UNFACT

10   10   5.997241E+05   1.667433E-06   1.667433E-06   1.667433E-06
11   11   5.997145E+05   1.667460E-06   1.667460E-06   1.667460E-06
12   12   5.993929E+05   1.668355E-06   1.668355E-06   1.668355E-06
13   13   7.098908E+05   1.667433E-06   1.667433E-06   1.667433E-06
14   14   7.098772E+05   1.667460E-06   1.667460E-06   1.667460E-06
...
4664   4664   1.377801E+06   1.389541E-06   1.389541E-06   1.389541E-06
4665   4665   1.376438E+06   1.390057E-06   1.390057E-06   1.390057E-06
4725   4725   0.000000E+00   0.000000E+00   0.000000E+00   0.000000E+00

for non-diagonal components:

                 FEAP_FACT   FEAP_UNFACT   PARD_FACT   PARD_UNFACT

Example of symmetric components:

11   14   3.939360E-01   6.568726E-07   6.568726E-07   6.568726E-07
14   11   3.939360E-01   6.568726E-07   6.568726E-07   6.568726E-07

Examples of non-symmetric components:

3836   3985   -8.093324E-08   -1.317032E-13   -1.317032E-13   -1.317032E-13
3985   3836   1.204221E-09   0.000000E+00   -1.317032E-13   -1.317032E-13
...
3910   3986   -5.151798E-10   0.000000E+00   -2.712450E-13   -2.712450E-13
3986   3910   -1.045583E-07   -2.712450E-13   -2.712450E-13   -2.712450E-13

My observations are:

1) In all the cases, there is no difference between Pardiso Factored and Unfactored tangents.

I enabled the debugger, Here are the outputs:
Pardiso -> UTANG

>>>>> PARDISO for non-symmetric indefinite matrices activated.

          Compressed Storage =   335223                    t=     0.18     0.00
 Reordering completed ...
Number of nonzeros in factors =            1
Number of factorization MFLOPS =            1
Factorization completed ...
Number of zero pivots        4311

Pardiso -> UTANG,,-1

>>>>> PARDISO for non-symmetric indefinite matrices activated.

          Compressed Storage =   335223                    t=     0.16     0.00
 *WARNING* Unfactored tangent produced do not try normal solution.

Everything looks fine in the output, so I think there is a problem in Factorization step of Pardiso.

2) in all the cases, tangent components of Pardiso are same as unfactored components of FEAP standard solver (just as I was suspecting).

3) obviously, UTANG command in Pardiso is not creating an unsymmetric tangent matrix (kindly have a look at the examples of non-symmetric components)

4) UTANG,,-1 just calculates the lower triangle of the non-symmetric tangent matrix.

Then, I tried to find if it is a problem with NURBFEAP or it happens even by FEAP itself.

I tried standard elements (and material models) by FEAP (attached input file) and again printed out the tangent matrices. I can confirm that the problem is with both FEAP and NURBFEAP. In all the cases, Pardiso calculates an unfactored symmetric tangent. The solution converges for simple problems (although the tangent is not exact, it can still lead the newton method to the correct solution. It is just the final residuals in each time step (although converged) are a little larger than the standard feap solver), but for my problem which is not symmetric, it produces a wrong tangent.

I hope this information helps to find the problem. (My element needs many modifications in feap (umacro, header files, user arrays, ...), I can send it to you, if necessary)

P.S: another thing which is interesting is that the size of output tangent are different, in the sense that there are components in UTANG (standard) that does not have any counterpart in Pardiso tangent. Kindly have a look at the last rows of attached data.rar->compare.txt

[Prof. R.L. Taylor]

Pardiso probably uses its own memory to do the factorization.  Are the answers the same?

[resammc]

It is actually very difficult to answer this question. Yes, the final answers seem to be the same. But, I did try something which I think should not affect the final result using different solvers:
I have this simple thermal element from FEAP Theory Manual modified to work with NURBFEAP (using nurbs shape functions,...). When I solve the 2D problem using 3 quadrature points in each direction (for a quadratic NURBS), the final results from FEAP Standard solver and Paradiso are same. But, when I reduce the quadrature points to just one in each direction, the final results are totally different (eg. 2.2106 in Paradiso vs. 4.7510 in FEAP solver).
So, I think the results are different, but the problem is not that complex to show this difference.

[Prof. R.L. Taylor]

But when you reduce to 1 quadrature point the tangent matrix is probably singular and thus the solution will be corrupted by round-off.    Doesn't the profile solver tell you at least 15 digits lost during factoring?

[resammc]

ah yes, I didn't check that, it does say that at least 7 digits lost in reducing diagonals of 1 equation and the maximum number of diagonal digits lost is 15.

[resammc]

Just to put a conclusion on this:

Special care must be taken in regards to the singularity of the tangent matrix when using Pardiso solver. In my problem (a solidification problem), since all the calculations are performed in the interphase between the liquid and solid phases, most of the nodes (control points) which are deep inside the liquid (or solid) phase will not participate in the calculations, hence large number of zero values in the tangent matrix diagonal. Profile solver is more flexible in handling these kind of tangents. For my case, I simply added some fake residual + tangent components for the not-participating degrees of freedom (+a corrector with history variables after each time step) and overcame the singularity problem. Now, Pardiso works as expected.

[Brisk]

I've encountered the same problem in FEAP 8.5 with pardiso solver. Different from resammc, the solution converges in the first server time steps, and the results is identical to FEAP's standard solver. However, after server steps, the program stopped with "Residual norm is NaN or Inf". If I use the stand solver, this error doesn't happen during the whole process.

When I used the standard solver, the output file tells me that "Maximum no. diagonal digits lost:  1". Is this the reason that pardiso solver doesn't converge?

OS: CentOS 6.5
Compiler+Math Library: Intel® Parallel Studio XE Cluster Edition for Linux* 2013