Problems with domain integration [v8.5]

[a118145]

Dear feap forum,

in one of my post processing routines, I am using isw=16 (JINT) to calculate something similar to a generalized J-integral. I am triggering the calculation using

Code: [Select]

call formfe(np(40),np(26),np(26),np(26),.false.,.true.,.false.,.true.,16,1,numel,1)in order to access the assembled FE-vector via hr(np(26)) in a subsequent routine. The results were qualitatively as expected, however the numbers, that I obtained using parfeap, were to low. This is where I got suspicious.

In order to debug the "almost correct" values, I changed the assembly in isw=16 in my user element in a way, that only the shapes themselves are numerically integrated and the results are stored in the first components of p, i.e.

Code: [Select]

p(1,1:nel) = p(1,1:nel) + Nscalar*w*detJNote, that the loop over the QPs and all other code lines are not displayed.

In a serial execution, the sum over all first component p-entries yields the surface area of my simulation domain, which is 1 in this case. Everything works as expected. However, using two (or more) partitions in parfeap, I obtain a reduced surface area of ~0.96 (which seems to be unaffected by the number of partitions). Note, that I am only accounting for the nodal contributions 1:numpn in case of parfeap. This is definitely not, what I expected. But what is even worse is the fact, that the sum over ALL nodes (including ghost nodes) 1:numnp in all partitions yields the same accumulated value of 0.96 as if it did not make any difference. To rule out any possible error which could be introduced by my MPI data collection routine, I even summed the individual contributions manually but this did not change anything.

I have no clue what is happening there but I suspect that something goes wrong during assembly. Does anybody have an idea?

Kind regards,

a118145

[Prof. R.L. Taylor]

It appears that ghost elements are only excluded when isw.eq.13, thus, I expect the problem is that the double counting may polute your answer.

[Prof. R.L. Taylor]

Is there an input file you can share -- and the test code you are using.

[a118145]

Dear Professor Taylor,

Thank you for your quick answer!

It appears that ghost elements are only excluded when isw.eq.13, thus, I expect the problem is that the double counting may polute your answer.

This would explain, why I see no difference between the summation over 1:numpn and 1:numnp. However, for this very simple case, I would expect a sum greater than 1 if there is a double counting of integrated surface areas.

Is there an input file you can share -- and the test code you are using.

I will try to set up a very simple test case. But before that, I will try to switch to isw = 13 in my user element and test the surface area calculation there. If that does the trick, we can confine the problem.

Kind regards,
a118145

[a118145]

Dear Professor Taylor,

I did not manage to trigger the calculation by exploiting isw = 13. Still, I am unsure about the "ghost element" statement you made. As far as I understand, the routine which invokes isw = 16 in my user element expects a nodal vector p, which is then assembled in the global vector. As there are only node related values in this vector, it should not be a question of active or inactive elements, but rather a matter of "ghost" or "not ghost" nodes. Up to now I thought, I have to make the latter decision in my personal accumulation routine, i.e. not counting ghost node contributions, because every ghost node has a corresponding real counterpart in the other partition.

I will try to supply you with a minimal working example including, which is not so easy because it is not just one file, but many, which rely on each other.

Kind regards,

a118145

[Prof. R.L. Taylor]

Sorry, it is not in yet.  There was a remark about this previously and it was noted that one had to skip the ghost elements in some instances.  I believe the correction was as indicated below in the routine pform.f in the parfeap directory:


!     Loop over active elements
      do n = nn1,nn2,nn3

        n_el = n

!       Check for active regions
        if((nreg.lt.0 .and. ix(nen1-1,n).ge.0)
     &                .or. (abs(ix(nen1-1,n)).eq.nreg)) then

!        Check to skip ghost elements
         if(ghostfl .and. ix(nen1-1,n).eq.21) cycle

The issue remains to defined the logical variable 'ghostfl' in an appropriate way for your use --  one way is to add before the loop
      if(isw.eq.13) then
        ghostfl = .true.
      else
         ghostfl = .false.
      endif
or whatever isw you want to use.

[a118145]

Dear Professor Taylor,

I set up a very simple mesh with three unit square elements side by side along the x-direction and played around a bit with integration of the surface area.

Serial run:
In the isw=16 task of my user element, I wrote the nodal coordinates and the corresponding contribution to the terminal. All results are as expected: The nodal values are equal to a quarter of the surface area for each element, which can be seen from the attached debug.txt, block starting at line 5.
Now, after triggering the assembly into p, cf. first post, by writing

Code: [Select]

call formfe(np(40),np(26),np(26),np(26),.false.,.true.,.false.,.true.,16,1,numel,1)I am calling another subroutine

Code: [Select]

call calcsurf(hr(np(26))which takes the assembled residual/reactions (pointer 26) and performs a loop over all nodes (no ghost nodes for serial run). Here, something strange happens. You find the output for increment 9 (of 10) in debug.txt, block starting at line 17. The first nodal contributions perfectly sum up to 3, which is the total surface area (correct), but this seems more or less a coincidence. I apply a compressive displacement along the x-direction, which yields some reaction forces. Now, these forces seem to be persistent in hr(np(26)). Somehow I expected, that this array would be zeroed at the very beginning of formfe(...). Hence, I added pzero before formfe, so that it reads

Code: [Select]

call pzero(hr(np(26)),ndf*numnp*2)
call formfe(...)
call calcsurf(...)
This did the trick for serial feap and removed some perturbations in the y-component (I did not mention them before in order to not distract). However, having added the pzero call, I now get a message at the very end of the successful simulation: Header field was corrupted. Ignoring this for a moment, I now switch to parfeap.

Parallel run:
At first, I implemented your suggestion with ghostfl in pform.f. Looking at the output within the isw = 16 task in my user element, I observe either 2 (flag set as you suggested) or 3 (permanently set to .false.) element calls with the correct nodal contributions (0.25 each). I expected this, but this seems not to solve my problem.

For better understanding, I attached a sketch of the partitioning. Circles are "real" nodes and crosses are "ghost" nodes. Switching to the output in calcsurf(...) again, I observe strange behavior. In partition 1 (cf. debug.txt line 54ff), the leftmost two nodes seem to be ignored by the assembly in formfe, because they are 0 instead of 0.25. The same can be seen in partition 2 (cf. debug.txt line 71ff) from nodes 3 and 4. All ghost node contributions are zero, which I also find weird, but I would have excluded them from summation anyways. Now, the lack of 4*0.25 of surface area explains, why I originally observed slightly smaller values in my original post. Some nodes are simply ignored during assembly.

Also for parfeap, I now receive a message

->Header field was corrupted
->corrupted size vs. prev_size
->Program received signal SIGABRT: Process abort signal.

at the very end of the simulation, but only with the pzero modification. Maybe, I made something wrong in call pzero(hr(np(26)),ndf*numnp*2).

Let me sum up some major questions/observations:
- The problem lies not within my isw = 16 task in the user element. Everything is calculated as expected and manual summation yields the correct total surface area (=3).
- hr(np(26)) has to be zeroed manually (?)
- In parfeap, nodal contributions seem to be lost, if at all they are assembled by formfe.
- It seems as if (columns of) ghost nodes are never assembled (at least in my case).

I will now try to put together a small user macro, that I can upload here, so that you can test it yourself.

Kind regards,

a118145

[a118145]

Dear Professor Taylor,

Please find attached a zip archive containing the necessary files for the MWE. In order to test parfeap, you have to modify the file I2d3elem as described therein. Furthermore, you have to add the supplied user macro and the subroutine calcsurf to your feap installation. The user macro is already called during solution (solve.2d3elem) using AXOS.

I tried to make it as simple as possible without any other routines, that I wrote. This, however, requires a few modifications in

- elements/solid2d/solid2d.f and
- elements/solid2d/sld2d1.f

Lines to be added/modified are marked with -->

solid2d.f
1) At around line 308, add the following:

Code: [Select]

!     Compute stress-divergence vector (p) and stiffness matrix (s)

      if(mech) then

-->     if(isw.eq.16) then
-->       !write(*,*) 'isw=16 called successfully'
-->       call sld2d1(d,ul,xl,ix,th,s,p,ndf,ndm,nst,isw,mther)
-->     endif
sld2d1.f
2) In line 107, add isw = 16 to the elseif

Code: [Select]

      elseif(isw.eq. 3 .or. isw.eq. 5 .or. isw.eq. 6 .or.
-->  &       isw.eq.13 .or. isw.eq.14 .or. isw.eq.16) then
3) Before the looping over the QPs starts, the p has to be zeroed. The formfe in the provided user macro seems to call every element twice. Frankly, I don't know why. It does not, when my user element is called.

Code: [Select]

!       Numerical integration loop
-->     if(isw.eq.16) then
-->       p(1:ndm,1:nel)=0.0d0
-->     endif
        nn = 0
        do l = 1,lint
4) Add the summation inside the numerical integration loop. I put it right before if(isw.eq.3 .or. isw.eq.6).

Code: [Select]

-->       if(isw.eq.16) then
-->         do k = 1,nel
-->           p(1,k) = p(1,k) + shp2(3,k,l)*jac(l)
-->         end do
-->         !if(l.eq.lint) then
-->         !  write(*,*) 'nodal: ', p(1,1:4)
-->         !endif
-->       endif

          if(isw.eq.3 .or. isw.eq.6) then
The commented sections are for easier debugging. Using the provided files and making the above modifications, I observe the same phenomena as in my previous post.

Thank you very much for your time, effort and consideration!

Kind regards,

a118145

[Prof. R.L. Taylor]

If you follow the logic in solid2d.f where you put the added statements you will see that after it comes back it will go right back into the element with calls below, hence the double call.  This is why you need to zero your p array as the second time it comes in it will not be zero.

It seems dangerous to interfere with an active isw in feap.  Why did you use 16 --use something like 60 to be sure you do not have unwanted calculations. 

If you give me a user element I will try it, but I do not want to edit existing elements, or use 16.

[a118145]

Dear Professor Taylor,

Please find attached a slightly modified zip archive with an additional user element. I also changed the isw to 60 and modified all routines accordingly. You will recognize, that elmt30.f is basically a copy of solid2d, except that I hardwired the numerical integration for isw = 60 directly within this file instead of calling another subroutine (lines 317ff).

For serial feap, everything works fine. For parallel feap, I observe the same odd behavior, that some nodal contributions do not appear, even if the elements are assembled (no exclusion of ghost elements). Long story short: the problem is persistent and reproducable.

Kind regards and thank you again for your indefatigable help

a118145

[Prof. R.L. Taylor]

I have modified your code a bit to recover the global nodes.  I also traced a problem to ./parfeap/pform.f  where a small change needs to be made.  In the assembly part (around line 459) change to read:

              if(jsw.ne.14) then
                if(pfeap_on) then
                  if (.not.efl) then                                      DELETE "jsw.eq.6 .and." ON THIS LINE
                    call pparlo(ld,mr(np(245)),ix(1,n),
     &                           mr(np(31)+ndf*numnp),nst,nen)
                    nar = nst
                  elseif (jsw.ne.8 .and. np(245).ne.0) then

There are changes to calcsurf & umacr24  - and a return in elmt30 to preclude any further checks after the isw.eq.60 part. 

Hopefully this will fix your problem.  I did not test for all options of equation forms nor several different b.c. form.  Also I did not collect all the sums to one location but the output should only show the global node numbers and the  associated integral.

[a118145]

Dear Professor Taylor,

First tests are working fine and yield expected results - also in combination with my results collection routine. Thank you very much! I will run a few more tests and if anything should not work as expected, I will come back to this thread.

EDIT: I took the modifications to my original problem and everything seems fine.

Warm regards,

a118145