Welcome, Guest
Username: Password: Remember me
Forgot your password? Forgot your username? Create an account
Forum
System Modules
TELEMAC-3D
Calling different mesh zones to assign head losses
Flat Threaded

TOPIC: Calling different mesh zones to assign head losses

Calling different mesh zones to assign head losses 6 years 7 months ago #29613

  • SDAC
  • SDAC's Avatar
Another thing to mention is that I use this in parallel, so would the node IDs need to be converted into local IDs to see the effects of the drag force?
The administrator has disabled public write access.

Calling different mesh zones to assign head losses 6 years 7 months ago #29614

  • c.coulet
  • c.coulet's Avatar
  • OFFLINE
  • Moderator
  • Posts: 3722
  • Thank you received: 1031
No because you don't use a specific node number in your program
Christophe
The administrator has disabled public write access.

Calling different mesh zones to assign head losses 6 years 7 months ago #29615

  • c.coulet
  • c.coulet's Avatar
  • OFFLINE
  • Moderator
  • Posts: 3722
  • Thank you received: 1031
Even you use double precision data you can't be sure of the last digit value.
In the informatics algorithm using real comparison is "forbidden"...

The order of the loop is not important here between 2D nodes and planes.
But you should optimize your algorithm.
You could start with the loop on 2D nodes and then enter in the plane loop.
Here you could compute I3D and NORM and choose the right value of AREA according the value of PRIVE2.
You finish with the computation of S1U, S1V and S1W and it's done

No needs to pass in the 2D nodes loop 3 times

Regards
Christophe
The administrator has disabled public write access.

Calling different mesh zones to assign head losses 6 years 7 months ago #29643

  • SDAC
  • SDAC's Avatar
Thankyou for your help!

The main problem is that the model runs fine but no effects of the drag are shown. I asked about the loop order is because I think the drag is only being applied on the bottom. Whilst slight differences may come from DP conversion I'm not seeing any actual drag results.

I've optimised the code (still requires some further work). I removed other code in the Fortran file to check there wasn't any confounding or unnoticed issues from elsewhere the the results are the same. Telemac reads the private variables fine. I looked at other case (on the forum and in T2D's dragfo.f)and I don't see an issue with the code applying vertical head losses to the nodes, but the issue seems to be that the drag equations aren't applied correctly.

I've attached the code, including the application of other head losses to four layers at the bottom. I think that these may also have the same problem.
CALL OS( 'X=0      ',X=S1U)
	     CALL OS( 'X=0      ',X=S1V)
      S1U%TYPR='Q'
	 S1V%TYPR='Q'
      IF(NONHYD) THEN
        S0W%TYPR='0'
        S1W%TYPR='Q'
	  CALL OS( 'X=0      ',X=S1W)
	ENDIF

	  DO I=1,NPOIN2
	   IF (PRIVE2D%ADR(1)%P%R(I).GT.0.1D0)THEN
   			DO IPLAN = NPLAN-33,NPLAN-30
				I3D = I+NPOIN2*(IPLAN-1)
				 IF  (IPLAN==(NPLAN-33))THEN
					AREFU = 1.6546D0
					AREFV = 1.9889D0
					AREFW = 13.8340D0

				ELSE IF (IPLAN==(NPLAN-32))THEN
					AREFU = 12.0348D0
					AREFV = 15.0840D0
					AREFW = 77.5191D0

				ELSE IF  (IPLAN==NPLAN-31) THEN
				 	AREFU = 12.3361D0
					AREFV = 14.8194D0
					AREFW = 76.6161D0

				ELSE IF  (IPLAN==NPLAN-30) THEN
					AREFU = 2.1330D0
					AREFV = 2.4954D0
					AREFW = 15.3330D0
			ENDIF

	    NORM =SQRT(UN3%R(I3D)**2+VN3%R(I3D)**2+WN3%R(I3D)**2)
			S1U%R(I3D)=0.5D0*AREFU*CD*NORM
			S1V%R(I3D)=0.5D0*AREFV*CD*NORM
			S1W%R(I3D)=0.5D0*AREFW*CD*NORM
			END DO
			END IF
			END DO

	DO I=1,NPOIN2
	IF (PRIVE2D%ADR(2)%P%R(I).EQ.2D0)THEN
   			DO IPLAN = NPLAN-29,NPLAN
				I3D = I+NPOIN2*(IPLAN-1)
			     AREA = 24.3816D0

	ELSE IF (PRIVE2D%ADR(2)%P%R(I).EQ.3D0)THEN
   			DO IPLAN = NPLAN-29,NPLAN
				I3D = I+NPOIN2*(IPLAN-1)
			     AREA = 77.6668D0

	
	ELSE IF (PRIVE2D%ADR(2)%P%R(I).EQ.4D0)THEN
   			DO IPLAN = NPLAN-29,NPLAN
				I3D = I+NPOIN2*(IPLAN-1)
			     AREA = 27D0

	 NORM =SQRT(UN3%R(I3D)**2+VN3%R(I3D)**2+WN3%R(I3D)**2)
			S1U%R(I3D)=0.5D0*AREA*CD2*NORM
			S1V%R(I3D)=0.5D0*AREA*CD2*NORM
			S1W%R(I3D)=0.5D0*AREA*CD2*NORM
			END DO
			END IF
			END DO
The administrator has disabled public write access.

Calling different mesh zones to assign head losses 6 years 7 months ago #29662

  • SDAC
  • SDAC's Avatar
I've noticed a couple of things:

1) Assigning head losses using private variables doesn't work - no head losses are assigned. Hard-coding for specific node IDs does. From the sortie file Telemac reads the private variables fine hwoever. I assume it's a coding issue. I'm using v7p2r3. Does this version require additional keywords/code for using private variables? Using the canal example and similar code to mine results in the images below:

With private variables:
novel.png

With hard-coded values:
vel.png


2) when head losses are applied successfully using the hardcoded version, I can only use the drag term for Ux and Uy - when I use the term for Uw I get the "segmentation fault" error. Why would this be? I've not managed to find out why Telemac is fine with parametrisation for Ux and Uy but not on Uw.

I used the code below:
DO I=1,NPOIN2
        IF(PRIVE2D%ADR(1)%P%R(I).GT.0) THEN ! using PRIVE
!        IF (I.EQ.555 .OR. I.EQ.551) THEN !hardcoded
          DO IPLAN=NPLAN-9,NPLAN-6 ! drag on the bottom
            I3D = I+NPOIN2*(IPLAN-1) ! 3D nodes
            NORM=SQRT(UN3%R(I3D)**2+VN3%R(I3D)**2+WN3%R(I3D)**2) 
			S1U%R(I3D)=0.5D0*AREAU*CD*NORM
			S1V%R(I3D)=0.5D0*AREAV*CD*NORM
!			S1W%R(I3D)=0.5D0*AREAW*CD*NORM
The administrator has disabled public write access.

Calling different mesh zones to assign head losses 6 years 7 months ago #29664

  • pilou1253
  • pilou1253's Avatar
  • OFFLINE
  • openTELEMAC Guru
  • Posts: 584
  • Thank you received: 106
Hi,

I have been away for some time.
It is strange that you can't manage to have it working because I use the same metholody often and it works fine.

You have copied several versions of your implementation that are slightly different so it's actually hard to follow how you really implemented your code. For example, in some parts you use PRIVE2D%ADR(1)%P%R(I) and in others PRIVE2D%ADR(2)%P%R(I).

Christophe pointed out that it is important to keep in mind that a test IF(A.EQ.1.D0) might not work since 1.D0 can be either 0.999...999 or 1.000...001. A solution is to use IF(ABS(A-1.D0).LT.1.D-6). Another thing that can maybe explain the issue is the test PRIVE2D%ADR(1)%P%R(I).GT.0 - it should be 0.D0 for double precision, I don't know if the compiler will automatically convert the integer 0 to a double precision.

Anyway, you are close to succeeding, it should only be a few things to correct.

PS: I realized that my comments in my earlier message were not exact, I misread your code. Sorry for that.

Good luck!

Best regards
PL
The administrator has disabled public write access.
The following user(s) said Thank You: SDAC

Calling different mesh zones to assign head losses 6 years 7 months ago #29669

  • SDAC
  • SDAC's Avatar
Hi,

Thanks for the advice Pilou.

I have used .GT.0D0 for implementing drag before to see if I could avoid the .EQ. issue and the head losses were still not applied. I've adopted the approach you suggested for the double precision problem and I've not had any luck.

Given that a) Telemac reads the variables fine, and b) the code for PRIVE matches the method others have used, perhaps it's something to do with the stored variables?
For the variables: I created a new Serafin variable (BOTTOM FRICTION); re-named the variable; mapped out polygon values onto the mesh. I'm not sure if there's another step I'm missing? I've attached the file in case anyone has time to look. The .f code is the less optimized version but it still runs fine when hard-coding points.

Additionally, is there an alternative to private variables that could work in paralleled? I saw in the Sysephe forum that FIND_VARIABLE could be used to read the variables from the geometry file. I assume this would be the appropiate thing to adapt? From there I could implement it in my code. For example, tentatively:
CALL FIND_VARIABLE('SERAFIN ',T2D_FILES(T2DGEO)%LU,'DRAG',
& ZR,MESH%NPOIN,IERR)
CALL OV('X=Y      ',ZR,DRAG%R,ZF,NPOIN)

DO I=1,NPOIN2
	    IF (ABS(ZR-4.D0).LT.4D-6)THEN

...

Many thanks!

File Attachment:

File Name: Roughness_test.zip
File Size: 162 KB
The administrator has disabled public write access.

Calling different mesh zones to assign head losses 6 years 7 months ago #29680

  • SDAC
  • SDAC's Avatar
Just realised the folder contained the incorrect .cas file...


File Attachment:

File Name: Test_case.cas
File Size: 4 KB
The administrator has disabled public write access.

Calling different mesh zones to assign head losses 6 years 7 months ago #29684

  • SDAC
  • SDAC's Avatar
Some progress. Using the code below head losses are assigned to the layers:
DO I=1,NPOIN2
	        DO IPLAN = NPLAN-13,NPLAN
				I3D = I+NPOIN2*(IPLAN-1)
			     
	   NORM =SQRT(UN3%R(I3D)**2+VN3%R(I3D)**2+WN3%R(I3D)**2)
			S1U%R(I3D)=0.5D0*PRIVE%ADR(2)%P%R(I3D)*CD2*NORM
			S1V%R(I3D)=0.5D0*PRIVE%ADR(2)%P%R(I3D)*CD2*NORM
			S1W%R(I3D)=0.5D0*PRIVE%ADR(2)%P%R(I3D)*CD2*NORM
			END DO
			END DO

However I get high velocities and elevation at the inlet and no flow everywhere else: the flow can't move down the channel. Nodes with no roughness elements have a value of 0, so I expected that there wouldn't be any drag applied at these nodes using the equation above. There appears to be enough resistance to prevent flow moving down the channel; drag is applied at every node.

I assumed that the drag equation would be multiplied by the values on I3D in the private variable, I guess this is wrong?


vel3.png
The administrator has disabled public write access.

Calling different mesh zones to assign head losses 6 years 7 months ago #29686

  • pilou1253
  • pilou1253's Avatar
  • OFFLINE
  • openTELEMAC Guru
  • Posts: 584
  • Thank you received: 106
Hi!

I had a quick look at your files but did not try to run the case.
I think your fortran file can be improved. You mix a lot of float/double precision/integers for what should be double precision variables. Once more, I don't know if this is the cause but it's not "clean".

I also noticed an error in your loop between NPLAN-13 and NPLAN-10: you have a test IF(NPLAN==(NPLAN-33)) - i assume you want 13 and not 33.

As Christophe also pointed out, you don't need to loop three times over NPOIN2, you can put all the three IF tests within the same loop.

Beware that PRIVE2D has dimension NPOIN2 (defined on the 2D mesh), so it's actually wrong to use the 3D node IDs in your latest code.

I paste below an example of how I add drag on nodes for which PRIVE2D(1) is greater than 0.1D0:
!*************PLLI drag for vegetation**************
      DO I=1,NPOIN2 
!       PRIVE2D%ADR(1)%P%R(I) = variable "VEGETATION" in geometry file defined as 1st 2D private array
        IF(PRIVE2D%ADR(1)%P%R(I).GT.0.1D0) THEN
          DO IPLAN=1,NPLAN 
            I3D = I+NPOIN2*(IPLAN-1) 
            NORM=SQRT(UN3%R(I3D)**2+VN3%R(I3D)**2+WN3%R(I3D)**2) 
            S1U%R(I3D)= 0.5D0 * N*DIAM/AREA * NORM * CD
            S1V%R(I3D)= 0.5D0 * N*DIAM/AREA * NORM * CD
!            S1W%R(I3D)= 0.5D0 * N*DIAM/AREA * NORM * CD
          ENDDO
        ENDIF
      ENDDO
!*************PLLI drag for vegetation**************

It works fine for me...

Good luck!

Best regards
PL
The administrator has disabled public write access.
Flat Threaded
Moderators: pham
Forum
System Modules
TELEMAC-3D
Calling different mesh zones to assign head losses

The open TELEMAC-MASCARET template for Joomla!2.5, the HTML 4 version.