Using "ENER" in parfeap vs in serial feap

[crepes]

Hi,

when using "ENER" after convergence, e.g. for a linear elastic 3D problem, then serial feap writes in the Ouput file something like

Code: [Select]

*Command  11 * ener                v:   0.00       0.00       0.00   
                                                           t=     3.57     0.00
 --> Energy for all materials

     Momentum and Energy for Problem
          Linear  Momentum   :  0.00000E+00  0.00000E+00  0.00000E+00
          Angular Momentum   :  0.00000E+00  0.00000E+00  0.00000E+00
          Kinetic   Energy   :  0.00000E+00
          Potential Energy   :  5.00000E+01
          Work by Loads      :  0.00000E+00
          Total     Energy   :  5.00000E+01
          Total Ang. Momentum:  0.00000E+00
 *Command  12 * next                v:   2.00       5.00       2.00   

while in parfeap it looks like

Code: [Select]

*Command  11 * ener                v:   0.00       0.00       0.00   
                                                           t=     1.77     0.40
 --> Energy for all materials
 *Command  12 * next                v:   2.00       5.00       2.00   

in each of (eight) output files

Should this behave like that? Or is my installation corrupt?

Kind regards

[Prof. S. Govindjee]

Looking at the code  (super quickly) it looks like the data exchanges to accumulate the results across processors has not been implemented but that does not explain why you not even seeing per-process outputs.  I'll have a deeper look but in the meantime perhaps you can add some debugging lines in program/pmacr1.f in the block of code that is supposed to write out the results from [ener].

[crepes]

Thank you for your reply.

Okay, I will try.

Maybe another question: In a user element, in serial feap, I was able to calculate those integrated values like the energy (or simply the volume) based on an allocated array in which I summed up the value of each element. I started at "n_el.eq.1" and stopped at "n_el.eq.numel" under isw.eq.13 where I printed the content of this allocated array entry into a file. And I allocated the array under isw.eq.1.

For parfeap, following the suggestion in http://feap.berkeley.edu/forum/index.php?topic=2343.0, I made a small change in parfeap/pform.f:

Code: [Select]

c         Check for ghost elements (no energy contribution)

          if (jsw.eq.13 .and. ix(nen1-1,n).eq.21) cycle

However, I still got a too big result (like 10% or so). Could the usage of numel as the maximum element number for summing up the quantities be a problem in parfeap?

[crepes]

Now not using the allocated array, but the "epl" array thing instead which is normally used for ENERg, and then using epl(8) to sum up the element volumes instead of normally the potential energy, I get the following behavior.


When I run a simulation with a cube having an edge length of 1 and having 16x16x16 elements, I get the following wrong value for the total volume:

Code: [Select]

1.0307617187500002


Dividing (1.0307617187500002 - 1.0) by the volume of one element, i.e. 1/(16x16x16)=0.000244140625 yields pretty much exactly 126.

So it includes some elements more than once I guess.

When I remove

Code: [Select]

c         Check for ghost elements (no energy contribution)

          if (jsw.eq.13 .and. ix(nen1-1,n).eq.21) cycle

from pform.f (cf. my last reply), I get

Code: [Select]

1.2478027343750002

which then yields 1015 elements when I divide it by the elemental volume.

To get these values, I just added

Code: [Select]

      write(*,*) epl(8)

to the code in pmacr1. Then I got 8 outputs (using 8 processes) which I summed up manually

[Prof. S. Govindjee]

I haven't had a chance to check yet but this indeed is the issue.   The partitioning is by node but to assemble the equations for a given node we include so called ghost elements.  If you look in the (parallel) input files the elements that you should ignore are marked to be in region 21.  Thus the summation loop when pfeap_on is true should skip elements where ix(nen1-1, n).eq.21, where n is the element number and nen1 is found  in sdata.h.

I hope to get to this next week along with some other bug fixes for parfeap that have been piling up but hopefully this is enough information for you to temporarily patch your version.

[Prof. S. Govindjee]

I just  read your post again, a bit more carefully.  It seems that the problem is more than just checking for region 21 flags.  It could be an initialization issue.

[crepes]

Thanks for your time. Unfortunately I had to decide between two plans:

1) Investigate feap routines to find out how to overcome this problem in an elegant way
2) Save all element values to a file and combine the element volumes using an external python script

I really needed to move on with my work, so I decided to approach plan "2"

If anyone is interested in my procedure for this:
-I deleted the modification in pform.f which I mentioned before
-using paraview, I noticed that some elements occur more than once in different partitions
-under isw=1, I created one file for each single process of parfeap
-under isw=13, I summed up the element volume using the gauss points, i.e. gauss quadrature
-also under isw=13, I calculated the centroid of the element using the XL array
-also under isw=13, I saved both, the element volume and the 3 coordinates of the centroid to the created file
-In a python script, I loaded all the files (in my case 8 because of 8 processes) and combined them which gave a huge matrix.
-each line represents one element and contains the centroid coordinates and the corresponding element volume
-I removed all duplicate lines using the numpy "unique" function.
-Summing up the remaining entries in the column of the element volumes resulted in the correct total volume

[Prof. S. Govindjee]

It is going to take me a few more days probably to get the next bug-fix release out.  In the meantime, you can try these files to see if they fix the problem for you. pform.f and pamcr7.F go in parfeap.  pmacr1.f goes in program.

[Prof. S. Govindjee]

Remark:  The routines that I posted change the way that ghost elements are determined and assigned when multiple processors need the same element for computing nodal equations.  The prior code was 99% correct but failed to flag some elements as ghost when they should have been.

The error did not affect any residual computations, but if you were using features like ENERgy to compute integrals over element domains for non-nodal quantities, then the numbers will have in most cases been slightly off due to some elements being double counted.

[crepes]

Thank you very much!

After two test simulations, it seems to have fixed my problem.

(Instead of replacing the files, I carefully made some changes line by line, since I am using version 85. However, I also have the code for 86 which made it easy to highlight the changes in your uploaded files using a software tool.)