Problem with cs_user_extra_operations

[frpageot]

Hello,

I am facing a strange (for me) problem that may have an explanation. I am rying to simulate a filter with a free volume before and after the inlet to prevent any interference of the results with those boundaries. I declared 2 internal faces as inlet and outlet of the filter and would like to calculate the pressure drop thru between them. I used then the cs_user_extra_operations to calculate the total pressure (including velocity) for each inegrated on the total surface. I did it previously on CS2 and it works pretty well with the usproj subroutine. I copy/paste and modify the code for CS3 and the new subroutine to have the same thing :

Code: Select all

allocate(lstelt(max(ncel,nfac,nfabor)))

Pinlet=0.D0
Sinlet=0.D0
Poutlet=0.D0
Soutlet=0.D0
Pinlet_chicane=0.D0
Sinlet_chicane=0.D0
DPchicane=0.D0
Poutlet_chicane=0.D0
Soutlet_chicane=0.d0


 call getfac('inlet_chicane', nlelt, lstelt)
  !==========

  do ilelt = 1, nlelt

    ifac = lstelt(ilelt)
    ii = ifacel(1,ifac)
    jj = ifacel(2,ifac)
    Sface=sqrt(surfac(1,ifac)**2 + surfac(2,ifac)**2 + surfac(3,ifac)**2)
    

    Pinlet_chicane = Pinlet_chicane + (rtp(ii,ipr)& 
                              +propfa(ii,ipprof(irom)) &
                             *0.5d0*(rtp(ii,iu)**2 + rtp(ii,iv)**2 &
                             +rtp(ii,iw)**2 ))*Sface


 


    Sinlet_chicane=Sinlet_chicane+Sface

 
  enddo
  
Pinlet_chicane=Pinlet_chicane/Sinlet_chicane


call getfac('outlet_chicane', nlelt, lstelt)
  !==========

  do ilelt = 1, nlelt

    ifac = lstelt(ilelt)
    ii = ifacel(1,ifac)
    jj = ifacel(2,ifac)
    Sface=sqrt(surfac(1,ifac)**2 + surfac(2,ifac)**2 + surfac(3,ifac)**2)
    

    Poutlet_chicane = Poutlet_chicane+(rtp(ii,ipr)++propfa(ii,ipprof(irom)) &
                             *0.5d0*(  rtp(ii,iu)**2 + rtp(ii,iv)**2 &
                            + rtp(ii,iw)**2 ))*Sface


    Soutlet_chicane=Soutlet_chicane+Sface

  enddo
  
Poutlet_chicane=Poutlet_chicane/Soutlet_chicane

DPchicane=Pinlet_chicane-Poutlet_chicane


write (nfecra, 1000) Pinlet, Pinlet_chicane, Poutlet_chicane, Poutlet, DPchicane
1000 format('         Pression entrée chicane = ', e14.5, /, &
            '         Pression sortie chicane = ', e14.5, /, &
            '         Delta P chicane= ', e14.5)

deallocate(lstelt)

and print the results in the listing. Once I launch calculation on 2 processes I have this error message :

MPI_ABORT was invoked on rank 0 in communicator MPI_COMM_WORLD
with errorcode 1.

NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
You may or may not see output from other processes, depending on
exactly when Open MPI kills them.
--------------------------------------------------------------------------
--------------------------------------------------------------------------
mpiexec has exited due to process rank 0 with PID 6740 on
node 0605-labolinu exiting improperly. There are two reasons this could occur:

1. this process did not call "init" before exiting, but others in
the job did. This can cause a job to hang indefinitely while it waits
for all processes to call "init". By rule, if one process calls "init",
then ALL processes must call "init" prior to termination.

2. this process called "init", but exited without calling "finalize".
By rule, all processes that call "init" MUST call "finalize" prior to
exiting or it will be considered an "abnormal termination"

This may have caused other processes in the application to be
terminated by signals sent by mpiexec (as reported here).
--------------------------------------------------------------------------
solver script exited with status 1.

Error running the calculation.

Check code_saturne log (listing) and error* files for details.


****************************
Saving calculation results
****************************

Error in calculation stage.

the error is then
Signal SIGFPE (exception en virgule flottante) intercepté !

Pile d'appels :
1: 0x401751 <cs_user_extra_operations_+0x6dd> (cs_solver)
2: 0x7f4ddb01950a <caltri_+0x34a6> (libsaturne.so.0)
3: 0x7f4ddafefa35 <cs_run+0xa35> (libsaturne.so.0)
4: 0x7f4ddafeef1a <main+0x14a> (libsaturne.so.0)
5: 0x7f4dda5ee76d <__libc_start_main+0xed> (libc.so.6)
6: 0x400fb9 <> (cs_solver)
Fin de la pile

I did try to launch the calculation only on one pocess and it worked without error but instead of having the result of the calculation I have the code "NaN". If I change the density to a value instead of a property (1.2 for instance) in this case I have the result of the calculation . If I launch the calculation without the cs_extra_operations sub, it works well with 2 processes. It seems that there is something wrong in my code but I am unable to find it.

[Yohann Eude]

Hello,

I just read your code and when I read the doc user_guide about irom:
irom: Property number corresponding to the density Stored at the cells and the boundary faces.

The density is available only in the cells and on boundary faces and not in internal faces, so try to change propfa(ii,ipprof(irom)) by propce(ii,ipproc(irom)).

Say me if it's Ok?

Yohann

[frpageot]

Hello,

Yes this is working now with 1 process, is I change my code as you wrote, Pressures are calculated. This is strange because I used like this with CS2.0 without any problem. Thanks for this.

But I still have the problem with 2 processes, still the same error message.

[Yvan Fournier]

Hello,

Your problem with 2 processes seems to be partition-dependent.

Don't forget to use parsom on the partial sums before dividing (check appropriate user examples if necessary).

If you did not use parsom with version 2.0, your code was incorrect. Either you were lucky (if all faces were on the same rank), or unlucky (if some faces were on one rank, others on the other, you did not have a division by 0, but your results were probably wrong).

Regards,

Yvan

[frpageot]

Hello Yvan,

I tried to use parsom as well, but it did not solve the problem. But I have to check witjh the irom on propc. With CS2, results were maybe wrong because I did not use parsom, but frankly they do not differ a lot with testings... I'll let you know, thanks.

[frpageot]

Hello,

OK, I found why I still had the error when using parsom. I was doing the parsom inside the do/enddo loop. Once done just before dividing, this is OK. Thanks.