open TELEMAC-MASCARET The mathematically superior suite of solvers

Hi Pham,

I am trying to understand the Fortran file of both porosity and drag force. In the file, which value is indication the porosity here below? Is it 0.2 value indication the porosity?
CALL OS( 'X=C ' , POROS , POROS , POROS , 1.D0 )
!
!

---

!
DO IELEM = 1 , NELEM
!
XX1 = ( X(IKLE%I(IELEM) )+
& X(IKLE%I(IELEM+NELMAX) )+
& X(IKLE%I(IELEM+2*NELMAX) ))/3.D0
YY1 = ( Y(IKLE%I(IELEM) )+
& Y(IKLE%I(IELEM+NELMAX) )+
& Y(IKLE%I(IELEM+2*NELMAX) ))/3.D0
!
IF(INPOLY(XX1,YY1,XSOM1,YSOM1,NSOM1)
& .OR. INPOLY(XX1,YY1,XSOM2,YSOM2,NSOM2)) THEN
POROS%R(IELEM) = 0.2D0 * ( 1.D0 + ABS(XX1/50.D0) )
ENDIF
!
ENDDO
!

Again, in the drag force fortran file, I am not getting the drag force and porosity value here. If you look at below, what is Por= 19.D0/20.D0 means? if it is porosity then where is drag force value here? Is it 19/20= 0.95 porosity?

POR=19.D0/20.D0
!
DO IELEM = 1 , NELEM
!
XX1 = ( X(IKLE%I(IELEM) )+
& X(IKLE%I(IELEM+NELMAX) )+
& X(IKLE%I(IELEM+2*NELMAX) ))/3.D0
YY1 = ( Y(IKLE%I(IELEM) )+
& Y(IKLE%I(IELEM+NELMAX) )+
& Y(IKLE%I(IELEM+2*NELMAX) ))/3.D0
!
IF(INPOLY(XX1,YY1,XSOM1,YSOM1,NSOM1)
& .OR. INPOLY(XX1,YY1,XSOM2,YSOM2,NSOM2)) THEN
IF(XX1.GE.-10.D0.AND.XX1.LE.10.D0) THEN
POROS%R(IELEM) = POR
ELSEIF(XX1.LT.-10.D0) THEN
POROS%R(IELEM) = POR - (1.D0-POR) * (XX1+10.D0) / 10.D0
ELSEIF(XX1.GT.10.D0) THEN
POROS%R(IELEM) = POR + (1.D0-POR) * (XX1-10.D0) / 10.D0
ENDIF
ENDIF
!
ENDDO ! IELEM

Hello,

I have not understood your question "which value is indication the porosity here below?".

Anyway, if you want to modify Fortran files, I strongly recommend you to learn Fortran and read the developer guide available here after:
wiki.opentelemac.org/doku.php?id=documentation_latest
Not too much time to explain line by line what is implemented in simple Fortan files.

Porosity field is stored in local variable POROS (argument of this subroutine). It is defined here depending on the abscissa of the element (here to be accurate, the center of gravity of the triangle):
- 19/20+1/20*(x-10)/10 or 19/20+1/20*(x+10)/10 or 19/20 and if inside a polygon, POR is a local variable of the subroutine.
- 0.2*(1+abs(x/50)) or 1 and if inside a polygon.

The drag force is not defined in the Fortran file dedicated to porosity (user_corpor) but in user_dragfo, stored in FUDRAG and FVDRAG that will be used for the source terms of equations.

Chi-Tuan

Thanks you for your reply. I am using the fortran subroutine for dragforce and I am using just one polygon point in the fortran file.
Each time I am getting the same error: "GRAGFO: Area of zone equal to zero"
Any idea why it can give this error? I am not sure where to work to fix this issue.

ITERATION 0 TIME: 0.0000 S
TELEMAC2D INITIALIZED
USING STREAMLINE VERSION 7.3 FOR CHARACTERISTICS
DRAGFO: AREA OF ZONE EQUAL TO ZERO
STOP 1

PLANTE: PROGRAM STOPPED AFTER AN ERROR
RETURNING EXIT CODE: 2

are you sure there is at least one element of your mesh fully inside the polygon?

I am sure

could you post a screenshot or give the files?

Hi, I added the fortran file here.

Ok

You define a new polygon but the old condition remains in the poros modification...

Here I attached the BOTTOM. Let me know if you need any other files