Divergence of thermal problem while having more DOFs

[JStorm]

I want to use parFEAP with an user element.
The element works fine in serial FEAP but not in parFEAP.
Then I have simplified the problem to a pure thermal problem while keeping the other DOFs.
Thus, I have a 2D element with 4 DOFs, where DOF 3 is the thermal problem.
The element tangent for DOFs 1, 2 and 4 are set to identity and the residuals to zero:

Code: [Select]

| I 0      0 |             | 0        |
| 0 K_temp 0 | * delta_u = | res_temp |
| 0 0      I |             | 0        |

This minimum working example is working in serial FEAP and in parFEAP with one CPU, but not in parFEAP with CPU > 1.

Could this be a bug in parFEAP or is there something wrong with my code?
I am using FEAP 8.6.1j.

[Prof. S. Govindjee]

This should work just fine.   

Also I do not see anything obviously wrong in the inputs or outputs.  (Note I did not examine the element).

The residual is the same in all 3 of your test cases at the first solve, so it looks like your residual computation is correct.
I would suggest to start by looking at the tangent matrix to make sure it is the same between the 1 and 2 process runs.

Also what solver are you using?  Perhaps you can try some other solvers? SuperLU, a simple CG?

[Prof. R.L. Taylor]

I have one suggestion for your element.  I would add a "save" to each subroutine to ensure that data is not lost. 

[Prof. R.L. Taylor]

Note also the "dinput.f" has been replaced by the function module "pinput.f" -- in fact dinput just calls pinput.

use:  var_inp = tinput(d(...), nn)

logicals

[Prof. R.L. Taylor]

I ran your problem using the feap element and get very different answers for the first time step.  I attach my file, perhaps there is an error in my set up.  Check carefully.

[Prof. R.L. Taylor]

Sorry, I had some bad b.c. and missed the density you used.  I get your answers with the serial version.  So it seems your tangent/residual are correct.

There are some rather dangerous Fortran coding issues where you use parameters to dimension arrays that are not formally allocated anywhere.  While this often works, it is not strictly correct.  You might try using "allocate" for these variables to help trace the difficulty in the parallel solution.

[Prof. R.L. Taylor]

I am looking at the input files for the parallel runs and note a "bug".  If you look at Idisc_000* for MATE you will find a Parameter  with no space between  the last material.

This is also in other files.  Can you try adding the blank line to check if that fixes the problem?

[Prof. R.L. Taylor]

Bug fix:  In ./program/pconst.f

In output to 'iwd'  add a "/" before the "a"  so line reads

      if(lopn) write(iwd,'(/a)') 'PARAmeters'

This will add a line to the output files for safety.

[Prof. R.L. Taylor]

That is not the bug -- you are missing the PETSc on command.

[JStorm]

Dear Prof. Taylor,

Thank you so much for spending so much time for investigating the problem.
I have also recognised that there was a missing empty line after the material block but I was sure that the correct material input is not disturbed by that and that the parameter line was not important.
Thank you for the fix.

I have tried it with PETSc,ON and the simulation is converging.
I will check the results on Monday and try to understand why this causes the issue.
So far, I thought that PETSc,ON is the default setting and I have run many simulation without this line.

Thank you for the programming hints, too.
I will carefully check the variable declarations in the original element code.

[JStorm]

There are some rather dangerous Fortran coding issues where you use parameters to dimension arrays that are not formally allocated anywhere.  While this often works, it is not strictly correct.  You might try using "allocate" for these variables to help trace the difficulty in the parallel solution.

Dear Prof. Taylor,

The complete user element is also working now when I use "PETSc,ON".

Please, can you give me a hint on the misstake with wrong array dimensions?
I only see that we have declared s(nst,nst) and p(nst) even though the dimensions differ in case of isw=8.

[Prof. R.L. Taylor]

In the subroutiine elementexample you dimension shp(nen)  and nr(ndm+1,nen) -- this is not strictly leagal, in my experience is works for small arrays but at some size may get lost or overwritten -- depending on compilers, computers, etc.

One can "allocate" them once using if( .not.allocated( shp) then) allocate shp(nen)

Just be cautious in the use --I am sure somewhere in the feap source you will find we still do it!!!

[Prof. R.L. Taylor]

One other comment, you should be able to put the allocate in the isw = 14 part -- it gets called only once for every element of the type, whereas isw = 3 gets called at every solution step.

[Prof. R.L. Taylor]

In looking at changes to make "petsc on" the default, I note that in ./parfeap/pdomain.F it is set.  However, it is not valid for all cases.  Can you try making the following changes if it will work for your problem when you omit the 'petsc on' command?
1. Add set: pfeap_dstr = .false. in pdomain.F

2. In pcontr where the call to pidreset add to read:

!       Form when all equations are included
        if(pfeap_blk .or. pfeap_bcin) then

You must have been quite fortunate that your previous problem solutions worked without the command!!!

[JStorm]

Dear Prof. Taylor,

The line "pfeap_dstr = .false." was already in pdomain.f.
I have add "pfeap_bcin" to pcontr.f and rerun the simulation without PETSc,ON, which works fine.

I do not know why all the previous simultions have worked.
The convergence behaviour was always as expected.
There is only one test case which has never worked and which I will try again.

Regarding the memory allocation based on subroutine arguments, I was not aware that this coding pattern is not covered by Fortran standard.
I have seen it so many times.

Thank you