Internal faces and boundary values

[Stéphane Blanchet]

Dear support and users,
 
I am facing a problem concerning the post-processing. Indeed, I would like to extract information about boundary values (for example the velocity at the outlet, etc.) and if possible use these information in my post-processing software (Paraview).
I also would like to extract data from internal faces. I used the "Interior faces selection for post-processing" and I expected additional post-processing files but it seems that only the listpre file contains more information.
 
Any idea would be very helpful.
 
Best Regards,
 
Stéphane

[David Monfort]

Hi Stéphane,
The option named "Interior face selection for post-processing" is only present for the preprocessor stage, and not for the solver stage. That's why you don't have the wanted behavior.
If you want to define specific zones for post-processing, then you will have to use advanced user files namely:
- usdpst.f90: definition of user post-processing zones and writers (i.e. the EnSight case) - usvpst.f90: definition of user post-processing variables - usmpst.f90: modification of user post-processing zones (to follow some species concentration e.g.) - usnpst.f90: modification of the output frequency
The first two subroutines are the more important.
There are numerous examples in each of these files but if you have any questions or trouble, feel free to ask us!
David

[Stéphane Blanchet]

Hi David,
Thank you for your answer. Now, I have no problem for selecting boundary faces. But  regarding the internal faces, I still have some problems. Indeed, it seems that Code_Saturne does not extract the internal faces since nothing appears in my post-processing software.  There is no problem with the color since I can see some information about those internal faces in the listpre file. So, I am guessing there is a problem with the code but I don't know where...
Here is my code for the extraction of the internal faces:

if (ipart .eq. 1) then     nommai = 'internal faces'           do ifac = 1, nfac       ! Adjacent cells       iel1 = ifacel(1,ifac)       iel2 = ifacel(2,ifac)       ! Adjacent cell colors (with only 1 property per cell)       icoul1 = iprfml(ifmcel(iel1),1)       icoul2 = iprfml(ifmcel(iel2),1)       ! Should the face belong to the extracted mesh ?       if ((icoul1.eq.312.and.icoul2.eq.313).or.(icoul1.eq.313.and.icoul2.eq.312)) then         nlfac = nlfac+1         lstfac(nlfac)= ifac       endif     enddo   endif

 
Thank you!
 
Stéphane

[David Monfort]

I don't see any obvious error here...
What you can do to check if everything works as intended to print the number of faces after your loop on the internal faces (don't forget to take into account the parallelism if needed ;-)). Something like will be enough:
nlfacg = nlfac
if (irangp.ge.0) then
  call parcpt(nlfacg)
endif
write(nfecra,*) "Total number of selected internal faces", nlfacg
If you still haven't found a solution, send us the complete usdpst.f90 so that I can check it.

[Stéphane Blanchet]

Thank you for the tip!
Actually, I used the wrong code. Now it works better with call getfac(...), I can see on the listing the total number of faces extracted. But at the end of the calculation, Code_saturne returns this error:

SIGFPE signal (floating point exception) intercepted!

I think it comes from the part of the code where I assign optional categories to user meshes (-1 for volume mesh, -2 for boundary mesh), since none of these categories fits for internal faces.
So, what is the best solution to export the same standard output?
Regards,
Stéphane
 

[David Monfort]

If I understand correctly, you want to add the interior faces user-mesh to the standard EnSight case ? Then you should use the following code in usdpst to associate your mesh (interior faces) to the default writer.
ipart = 1   ! the mesh
icas = -1  ! the default writer
call pstass(ipart, icas)
Let me know if it doesn't work as intended.
Perhaps could you attach the usdpst and usvpst files so that I can check them.
David

[Stéphane Blanchet]

Well, I would like to export the internal face with the same output data as the standard EnSight case (velocity, pressure, etc.), but in a different file.

It works fine with the boundary faces with this code:

 ipart = 2  icat = -2  call pstcat(ipart, icat)

But it does not work with this code for the internal face:

 ipart = 1
 icat = -1
 call pstcat(ipart, icat)

Let me know if it is still not clear enough.

Here are my files (I did not change the usvpst file).

Stephane

[David Monfort]

You cannot do that easily, and not with the pstcat subroutine. With this subroutine, all you can do is associating variables defined on a given mesh with another mesh of the same kind of elements (cells -> cells, interior faces -> interior faces, boundary faces -> boundary faces).
Here, you want to associate variables defined on cells with a mesh containing only interior faces... that is not possible :) The solution is then to use the usvpst.f90 subroutine and mimick one of the different examples.

[Stéphane Blanchet]

Thank you, that's what I thought concerning the pstcat subroutine. Now I use the example 1.3 in the usvpst.f90 subroutine (pressure computed by interpolation, etc.) and it works fine.
But when I check the result on Paraview, and compare it with a "slice" of the main volume mesh (in order to check the same internal face), the pressure computed by interpolation and the pressure on the "slice" are quite different (about 10% of difference for the mean pressure).
I guess Paraview does not really take the internal faces values when doing a "slice", but 10% seems quite big to me...
One more question: In the usvpst.f90 subroutine, how can I synchronized neighbors for parallel calculations?
Regards,
Stéphane

[David Monfort]

I don't know about ParaView... but IMHO, the Code_Saturne calculation will always be better than the ParaView one because you can control the interpolation.
For parallel and periodic synchronisation, you should refer to the examples in usproj.f90.
As an example, here is the synchronisation for a scalar (keep the order of the subroutine calls):
! Parallel synchronization
if (irangp.ge.0) then
  call parcom(var)
endif

! Periodic synchronization
if (iperio.eq.1) then
  idimte = 0
itenso = 0
  call percom(idimte, itenso,
var, var, var,
              var, var, var,
              var, var, var)
endif
and a vector:
! Parallel synchronization
if (irangp.ge.0) then
  call parcom(varX)
  call parcom(varY)
  call parcom(varZ)


endif

! Periodic synchronization
if (iperio.eq.1) then
  idimte = 1
  itenso = 0
  call percom(idimte, itenso,
varX, varX, varX,
              varY, varY, varY,
              varZ, varZ, varZ)
endif