getfbr with Syrthes coupling

[FannyD]

Hello,

I’d like to know if there is a function working like ‘getfbr’ when coupling Code_Saturne with Syrthes.
In fact, I wrote a function computing the Nusselt number on the boundary between solid and fluid meshes. It works when starting Code_Saturne only. But when I launch a coupled calculation, the ‘getfbr’ can’t find any cell.

I could use ‘getcell’ with geometric conditions, but I’d like to have a code working with any mesh or any geometry (the height of my solid obstacle is variable).
Is there a way to use ‘getfbr’ here? Or another function working the same way?

Here is my code:

Code: Select all

integer iscal, ivar, irangv
integer ifac, iel, nlelt, ielcote, ilelt
integer, allocatable, dimension(:) :: lstelt 

double precision, allocatable, dimension(:,:) :: Nux
double precision Sloc, Tcell, Tcote, gradT, Temp0
double precision yy, zz, xx, delta, Nusselt, posi
integer nb_tri, ii, jj
integer impout1

!===============================================================================

!===============================================================================
! Initialization
!===============================================================================

delta = 10.d0/1000.d0
Temp0 = 0.d0
nb_tri = 0

allocate(lstelt(nfabor))  ! temporary array for boundary faces selection
!===============================================================================
! User operations
!===============================================================================

! On calcule le Nusselt autour de l'obstacle

if (ntcabs.eq.ntmabs) then  ! ne calcule qu'au dernier pas de temps

  iscal = iscalt         ! temperature scalar number
  ivar =  isca(iscal)    ! temperature variable number
  
  ! Appel des cellules de bord
  call getfbr('parois_obstacle',nlelt,lstelt) 
  allocate(Nux(nlelt,2))
  Nux = 0.d0
  
  do ilelt = 1, nlelt 
    ifac = lstelt(ilelt)
    iel = ifabor(ifac) ! element contre la cellule
    xx = xyzcen(1,iel)
    yy = xyzcen(2,iel)
    zz = xyzcen(3,iel) 

    ! Détermination de la position de la cellule
    if (xx.lt.2) then
      posi = yy
    else if (xx.gt.2.25) then
      posi = 0.5+(0.25-yy)
    else
      posi = 0.25+(xx-2)
    endif
    
    ! recherche de la cellule juste à côté de la cellule considérée.
    if (xx.lt.2) then
      call findpt(ncelet, ncel, xyzcen, xx-delta, yy, zz, ielcote, irangv)
    else if (xx.gt.2.25) then
      call findpt(ncelet, ncel, xyzcen, xx+delta, yy, zz, ielcote, irangv)
    else
      call findpt(ncelet, ncel, xyzcen, xx, yy+delta, zz, ielcote, irangv)
    endif
    

    !calcul du gradient de temperature
    Tcote = rtp(ielcote,ivar)
    Tcell = rtp(iel,ivar)  
    gradT = abs((Tcell-Tcote)/delta)
    ! calcul du nombre de Nusselt
    Nusselt = gradT/(Tcell-Temp0)
    !Nusselt = 1 !test pour vérifier que toutes les cellules sont visitées
    
    ! Mise a jour du vecteur du Nusselt : insertion de l'élément de manière triée
    if (nb_tri.eq.0) then 
      Nux(1,1) = posi
      Nux(1,2) = Nusselt
      nb_tri = nb_tri+1
    else
      ii = 1
      do while (posi.lt.Nux(ii,1))
        ii = ii+1
      enddo
      do jj = nb_tri+1,ii,-1
        Nux(jj,1) = Nux(jj-1,1)
        Nux(jj,2) = Nux(jj-1,2)
      enddo
      Nux(ii,1) = posi
      Nux(ii,2) = Nusselt
      nb_tri = nb_tri+1
    endif
 enddo
  

  ! Ecriture
  if (irangp.le.0) then
    open(impout1,file='Nusselt.dat', position = "append")
    write(impout1,*)  &
         '# abscisse_curviligne       Nusselt (x) '
    do ii=1,nlelt
      write(impout1,120)  &
        Nux(ii,1), Nux(ii,2)
    enddo
    close(impout1)
  endif

deallocate(lstelt)
deallocate(Nux)

endif

Thanks in advance

[Yvan Fournier]

Hello,

There is no reason getfbr should not work when coupling with Syrthes.

But using findpt in a loop on elements or faces is a bad idea, as it will cause a failure most of the time when running in parallel (it is a collective call, and must be called the same number of times by each rank).

There might be more generic solutions to your problem. Did you look at the various user examples or functions in post_util.f90 to compute a Nusselt number at the cells just adjacent to the boundary ?

Otherwise, I recommend adding logging (print) or better, running under a debugger to understand what is going wrong in the coupled case.

Regards,

Yvan

[FannyD]

Hello,

Concerning the calculation of the Nusselt number, I'm using using this formula to stay as close as possible of Results I get from an article. Once I get good enough results this way, I'll take a look at the the calculation of Nusselt using wall temperature or temperature gradient (on the basis of the cs_f_user_extra_operations subroutine of doxygen, as I could'nt find any file named post_util.f90, either on the SRC directory, or on the internet).

For the getfbr function, I first want to make it work on one proc (I've already an another function giving good results for the Nusselt number in parrallel, but without Syrthes coupling). For parrallel calculations, I effectively get failure when keeping findpt in a loop over elements, an I already correct this.
So as my sub work well (parrallel or not) without Syrthes, I tried to isolate my problem, running just the following code :

Code: Select all

integer nlelt

 impout1 = impusr(1)
  if (irangp.le.0) then
    open(impout1,file='test.dat', position = "append")
  endif

  call getfbr('parois_obstacle',nlelt,lstelt) 

  if (irangp.le.0) write(impout1,*)  nlelt

  if (irangp.le.0) close(impout1)

Without Syrthes, I get the correct number (42). With Syrthes, I get 0.
For the coupling, I followed the tutorial 'Three 2D disks', adapting it to my case. The results show a good agreement with my reference article, so the coupling seems to be correct.

I thought that maybe Code_saturne didn't see the surface between fluid and solid domains as a boundary? (getfbr function worked correctly when calling inlet/outlet surfaces).
Or could this come from a bad installation of Syrthes/a missing library?

I'm working with Code_Saturne 3.0.7 and Syrthes 4.1.1.

Regards,

Fanny

[Yvan Fournier]

Hello,

On how many parallel ranks are you using each code ?

The post_util.f90 file is in the code's sources (under src/base).

Also, Nusselt number examples may be found in SRC/EXAMPLES.

Regards,

Yvan

[FannyD]

Hello,

thanks for the location of post_util file.

For the test I made on getfbr, I'm running each code on one rank (one for C_S, one for Syrthes).
For the calculation, I run it with one rank for Syrthes and 15 for Code_Saturne.

Regards,

Fanny

[Yvan Fournier]

Hello,

Note that when running on more than 1 rank, the values returned by getfbr are local, but only rank 0 prints info in the listing by default. So to check your count, you should use something of the form:

integer nelt, neltg
...

Code: Select all

call getfbr('parois_obstacle',nlelt,lstelt) 
neltg = nelt
if (irangp.ge.0) call parcpt(neltg)

write(nfecra, *) 'selected faces: ', neltg

And you must avoid the use of indpt in loops on faces (oe other mesh-dependent elements)

Regards,

Yvan