|
Welcome,
Guest
|
|
TOPIC: Assigning vertical head losses
Assigning vertical head losses 7 years 11 months ago #27549
|
Hello!
I've written a code to model head losses at the first four bottom layers of a channel. It's heavily inspired by Pilou's approach to surface head losses(opentelemac.com/index.php/assistance/for...log-boom-in-3d#26925) My head loss equation has varying area values depending on U/V/W, and depending on which layer. I've written a code that identifies a layer and then defines the area values for U/V/W for that particular layer, before inputting them into my head loss equation. The idea is to have each layer to model a different head loss for U/V/W as suggested by the area. However, it's not working. I think this comes from my rusty knowledge of Fortran. Is the below appropriate for assigning different area values per U/V/W per particular layer? AREFU/AREFV/AREFW are the area values to be assigned for layers NPLAN-17, NPLAN-16, NPLAN-15, NPLAN-14 respectively. DO I=1,NPOIN2
IF (PRIVE%ADR(1)%P%R(I).GT.0.1D0)THEN !search for points
! DO IPLAN = NPLAN-17,NPLAN-14
I3D = I+NPOIN2*(IPLAN-1)
IF (NPLAN=(NPLAN-17))THEN
AREFU=1.6546
AREFV=1.9889
AREFW=13.8340
ELSE IF (NPLAN=(NPLAN-16))THEN
AREFU= 12.0348
AREFV=15.0840
AREFW=77.5191
ELSE IF (NPLAN==NPLAN-15) THEN
AREFU= 12.3361
AREFV=14.8194
AREFW=76.6161
ELSE IF (NPLAN==NPLAN-14) THEN
AREFU= 2.1330
AREFV=2.4954
AREFW=15.3330
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
ENDIF
END DOMany thanks! |
|
The administrator has disabled public write access.
|
Assigning vertical head losses 7 years 11 months ago #27551
|
Hello,
Your code seems good except for two details: First: On the lines 5 and 9 of the piece of code you posted (IF (NPLAN=(NPLAN-17))THEN and ELSE IF (NPLAN=(NPLAN-16))THEN) the conditional is incorrect. The = sign means an assignment, not a logical test. You have to replace by a ==, that will perform the conditional test. Second, after you fix the above, the conditionals say: IF (NPLAN==(NPLAN-17)). But NPLAN will never be equal NPLAN-17. Probably you mistyped the variable name here. What you probably want is IF (IPLAN==(NPLAN-17)). I believe that if you change these bits your code should run fine. Hope it helps. Best regards, Phelype |
|
The administrator has disabled public write access.
The following user(s) said Thank You: SDAC
|
Assigning vertical head losses 7 years 11 months ago #27552
|
Thanks for you help, Phelype!
However, I now get the error: "segmentation fault - invalid memory reference." I used the DEBUGGER = 1, and the backtrace gives: #0 ffffffffffffffff I've not been able to figure out where this is stemming from. I don't suppose anyone has ideas how to overcome this? I've not had a problem with this before, so I imagine it might be something to do with my fortran file? Loading Options and Configurations
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
_____ __ __
|___ | /_ | /_ |
__ __ _/ / _ __ | | _ __ | |
\ \ / / |_ _| | '_ \ | || '__| | |
\ V / / / | |_) | | || | | |
\_/ /_/ | .__/ |_||_| |_|
| |
|_|
_ _ __ ___ ___ _ _ ___
_| || |_ /_ | / _ \ |__ \ | || | / _ \
_ __ ___ __ __ |_ __ _| | || | | | ) || || |_ | | | |
| '__| / _ \\ \ / / _| || |_ | || | | | / / |__ _|| | | |
| | | __/ \ V / _ |_ __ _| | || |_| | / /_ | | | |_| |
|_| \___| \_/ (_) |_||_| |_| \___/ |____| |_| \___/
... parsing configuration file: C:\opentelemac-mascaret\v7p1\configs\systel_v7p1.cfg
Running your CAS file for:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+> configuration: win7gfors
+> root: C:\opentelemac-mascaret\v7p1
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
... reading the main module dictionary
... processing the main CAS file(s)
+> running in English
... handling temporary directories
... checking coupling between codes
... checking parallelisation
... first pass at copying all input files
copying: PHEOX_BC.cli C:\TELEMAC\PH04 Telemac\PHEOX_KS0.01PSI.cas_2017-08-14-19h16min48s\T3DCLI
copying: PHEOX.slf C:\TELEMAC\PH04 Telemac\PHEOX_KS0.01PSI.cas_2017-08-14-19h16min48s\T3DGEO
copying: Drag2.f C:\TELEMAC\PH04 Telemac\PHEOX_KS0.01PSI.cas_2017-08-14-19h16min48s\t3dfort.f
re-copying: C:\TELEMAC\PH04 Telemac\PHEOX_KS0.01PSI.cas_2017-08-14-19h16min48s\T3DCAS
copying: telemac3d.dico C:\TELEMAC\PH04 Telemac\PHEOX_KS0.01PSI.cas_2017-08-14-19h16min48s\T3DDICO
... checking the executable
re-copying: Drag2.exe C:\TELEMAC\PH04 Telemac\PHEOX_KS0.01PSI.cas_2017-08-14-19h16min48s\out_Drag2.exe
... handling sortie file(s)
Running your simulation(s) :
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C:\TELEMAC\PH04 Telemac\PHEOX_KS0.01PSI.cas_2017-08-14-19h16min48s\out_Drag2.exe
_ _ ____ _ _____ __
| | | | |___ \ | | |___ | /_ |
| |_ ___ | | ___ _ __ ___ __ _ ___ __) | __| | ______ __ __ _/ / _ __ | |
| __| / _ \| | / _ \| '_ ` _ \ / _` | / __| |__ < / _` | |______| \ \ / / |_ _| | '_ \ | |
| |_ | __/| || __/| | | | | || (_| || (__ ___) || (_| | \ V / / / | |_) | | |
\__| \___||_| \___||_| |_| |_| \__,_| \___||____/ \__,_| \_/ /_/ | .__/ |_|
| |
|_|
------------------------------------------------------------------------------
LISTING OF TELEMAC-3D
TTTTT EEEEE L EEEEE M M AAAAA CCCCC
T E L E MM MM A A C
T EEE L EEE M M M AAAAA C
T E L E M M A A C
T EEEEE LLLLL EEEEE M M A A CCCCC
3D VERSION 7.2 FORTRAN 90
********************************************
* LECDON: *
* BEFORE CALLING DAMOCLES *
********************************************
THE KEY-WORD: PARALLEL PROCESSORS
APPPEARS AT LEAST TWICE , THE LAST VALUE WILL BE KEPT...
END OF FILE FOR DAMOCLES
********************************************
* LECDON: *
* AFTER CALLING DAMOCLES *
* CHECKING OF DATA READ *
* IN THE STEERING FILE *
********************************************
EXITING LECDON. NAME OF THE STUDY:
TELEMAC 3D : Channel
OPENING FILES FOR TELEMAC3D
POINT_TELEMAC3D: MEMORY ALLOCATION
READ_MESH_INFO: TITLE= PheoX
NUMBER OF ELEMENTS: 31269
NUMBER OF POINTS: 16341
FORMAT NOT INDICATED IN TITLE
MXPTEL (BIEF) : MAXIMUM NUMBER OF ELEMENTS AROUND A POINT: 8
MAXIMUM NUMBER OF POINTS AROUND A POINT: 8
CORRXY (BIEF):NO MODIFICATION OF COORDINATES
MESH: MESH2D ALLOCATED
READ_MESH_INFO: TITLE= PheoX
NUMBER OF ELEMENTS: 31269
NUMBER OF POINTS: 16341
FORMAT NOT INDICATED IN TITLE
MXPTEL (BIEF) : MAXIMUM NUMBER OF ELEMENTS AROUND A POINT: 8
MAXIMUM NUMBER OF POINTS AROUND A POINT: 8
CORRXY (BIEF):NO MODIFICATION OF COORDINATES
MESH: MESH3D ALLOCATED
2D MESH
-------
2D ELEMENT TYPE : 10
NUMBER OF 2D NODES : 16341
NUMBER OF 2D ELEMENTS : 31269
NUMBER OF 2D BOUNDARY NODES : 1411
3D MESH
-------
3D ELEMENT TYPE : 40
NUMBER OF 3D NODES : 245115
NUMBER OF 3D ELEMENTS : 437766
NUMBER OF LEVELS : 15
NUMBER OF BOUNDARY ELEMENTS : 19754
TOTAL NUMBER OF BOUNDARY NODES : 53847
INCLUDING LATERAL BOUNDARIES : 21165
SURFACE : 16341
BOTTOM : 16341
END OF MEMORY ALLOCATION
APPEL DE LECLIM
RETOUR DE LECLIM
APPEL DE INBIEF POUR MESH2D
INBIEF (BIEF): NOT A VECTOR MACHINE (ACCORDING TO YOUR DATA)
RETOUR DE INBIEF
APPEL DE INBIEF POUR MESH3D
RETOUR DE INBIEF
CHECK: NO ERROR HAS BEEN DETECTED
APPEL DE FONSTR
VARIABLE DRAG FOUND IN THE GEOMETRY FILE
STRCHE (BIEF): NO MODIFICATION OF FRICTION
RETOUR DE FONSTR
CALLING T3D_CORFON
RETURN FROM T3D_CORFON
APPEL DE FRONT2
THERE IS 2 LIQUID BOUNDARIES:
BOUNDARY 1 :
BEGINS AT BOUNDARY POINT: 679 , WITH GLOBAL NUMBER: 16272
AND COORDINATES: 122.5674 28.19508
ENDS AT BOUNDARY POINT: 708 , WITH GLOBAL NUMBER: 16335
AND COORDINATES: 123.1460 33.40766
BOUNDARY 2 :
BEGINS AT BOUNDARY POINT: 1388 , WITH GLOBAL NUMBER: 25
AND COORDINATES: 12.79237 37.78787
ENDS AT BOUNDARY POINT: 1 , WITH GLOBAL NUMBER: 1
AND COORDINATES: 12.78001 33.01963
THERE IS 2 SOLID BOUNDARIES:
BOUNDARY 1 :
BEGINS AT BOUNDARY POINT: 1 , WITH GLOBAL NUMBER: 1
AND COORDINATES: 12.78001 33.01963
ENDS AT BOUNDARY POINT: 679 , WITH GLOBAL NUMBER: 16272
AND COORDINATES: 122.5674 28.19508
BOUNDARY 2 :
BEGINS AT BOUNDARY POINT: 708 , WITH GLOBAL NUMBER: 16335
AND COORDINATES: 123.1460 33.40766
ENDS AT BOUNDARY POINT: 1388 , WITH GLOBAL NUMBER: 25
AND COORDINATES: 12.79237 37.78787
RETOUR DE FRONT2
APPEL DE CONDIM
CONDIM: DYNAMIC PRESSURE INITIALISED TO ZERO
HYDROSTATIC PRESSURE INITIALISED TO ZERO.
RETOUR DE CONDIM
APPEL DE MESH_PROP
RETOUR DE MESH_PROP
APPEL DE VERMOY
RETOUR DE VERMOY
APPEL DE LICHEK
RETOUR DE LICHEK
APPEL DE MASBAS2D
RETOUR DE MASBAS2D
APPEL DE GRAD2D
RETOUR DE GRAD2D
APPEL DE FSGRAD
RETOUR DE FSGRAD
APPEL DE KEPINI
RETOUR DE KEPINI
APPEL DE DRSURR
RETOUR DE DRSURR
APPEL DE COEFRO
RETOUR DE COEFRO
APPEL DE TFOND
RETOUR DE TFOND
APPEL DE VISCKE
RETOUR DE VISCKE
================================================================================
ITERATION 0 TIME 0 D 0 H 0 MN 0.0000 S ( 0.0000 S)
================================================================================
APPEL DE PRERES_TELEMAC3D
RETOUR DE PRERES_TELEMAC3D
APPEL DE WRITE_HEADER
RETOUR DE WRITE_HEADER
APPEL DE WRITE_MESH
RETOUR DE WRITE_MESH
APPEL DE BIEF_DESIMP
RETOUR DE BIEF DESIMP
APPEL DE WRITE_HEADER
RETOUR DE WRITE_HEADER
APPEL DE WRITE_MESH
RETOUR DE WRITE_MESH
APPEL DE BIEF_DESIMP POUR 2D
RETOUR DE BIEF_DESIMP POUR 2D
--------------------------------------------------------------------------------
MASS BALANCE
INITIAL MASS OF WATER IN THE DOMAIN : 278.44066123497493
PREMIER APPEL DE SOURCES_SINKS
RETOUR DE SOURCES_SINKS
PREMIER APPEL DE PREADV
RETOUR DU PREMIER APPEL DE PREADV
BOUCLE EN TEMPS LT= 1
APPEL DE VERMOY
RETOUR DE VERMOY
APPEL DE COEFRO
RETOUR DE COEFRO
APPEL DE LICHEK
RETOUR DE LICHEK
APPEL DE BORD3D
RETOUR DE BORD3D
APPEL DE TBORD
RETOUR DE TBORD, APPEL DE TFOND
RETOUR DE TFOND
APPEL DE KEPCL3
RETOUR DE KEPCL3
APPEL DE TRISOU
RETOUR DE TRISOU, APPEL DE SOURCE
_____________
runcode::main:
:
|runCode: Fail to run
|C:\TELEMAC\PH04 Telemac\PHEOX_KS0.01PSI.cas_2017-08-14-19h16min48s\out_Drag2.exe
|~~~~~~~~~~~~~~~~~~
|Program received signal SIGSEGV: Segmentation fault - invalid memory reference.
|
|Backtrace for this error:
|#0 ffffffffffffffff
|~~~~~~~~~~~~~~~~~~ |
|
The administrator has disabled public write access.
|
Assigning vertical head losses 7 years 11 months ago #27554
|
The segmentation fault happens when you do something incorrect that even the compiler didn't expect.
First you need to find out exactly where this is happening (the line of code with a problem), then find out why it is happening. What I do when this kind of thing happens is put a numbered PRINT*, statement after every line to check which one is throwing the segmentation fault. Something like this: PRINT*, 'I AM HERE 1'
DO I=1,NPOIN2
PRINT*, 'I AM HERE 2'
IF (PRIVE%ADR(1)%P%R(I).GT.0.1D0)THEN !search for points
PRINT*, 'I AM HERE 3'
! DO IPLAN = NPLAN-17,NPLAN-14
PRINT*, 'I AM HERE 4'
I3D = I+NPOIN2*(IPLAN-1)
PRINT*, 'I AM HERE 5'
IF (NPLAN==NPLAN-17) THENAlso, I just noticed that in the piece of code you posted, the line DO IPLAN = NPLAN-17,NPLAN-14 is commented. If it is also commented in your code, then your error might be coming from the extra END DO. I think this line should not be commented... Best regards, Phelype |
|
The administrator has disabled public write access.
The following user(s) said Thank You: SDAC
|
Assigning vertical head losses 7 years 11 months ago #27557
|
Hi there!
I am in a rush and did not take the time to look at your code in depth but I have a quick remark: did you set NUMBER OF PRIVATE ARRAYS to 1 or larger in your steering file? If not then this could explain the segmentation fault since the structure PRIVE%ADR(1)%P%R(I) has not been allocated. Another linked remark: PRIVE%ADR(1)%P%R(I) refers to the 3D private arrays, if you do as my method, looking into a variable stocked in the geometry file, this is actually a 2D private variable and you need to use NUMBER OF 2D PRIVATE ARRAYS and PRIVE2D%ADR(1)%P%R(I). So the test: DO I=1,NPOIN2 IF(PRIVE%ADR(1)%P%R(I)... don't seem logical to me (but I can be wrong since I dont know what you actually put in PRIVE%ADR(1)%P%R(I)). Good luck, PL |
|
The administrator has disabled public write access.
The following user(s) said Thank You: SDAC
|
Assigning vertical head losses 7 years 11 months ago #27560
|
Hi,
a good way to detect segmentation faults is to use some more strict compiler options during the development. For gfortran you can use -fbounds-check -fbacktrace Good luck Leo |
|
The administrator has disabled public write access.
|
Assigning vertical head losses 7 years 10 months ago #27610
|
Thanks for the help everyone. Removing the commented section allowed the code to run. It still ran without changing the arrays from RIVE%ADR(1)%P%R(I) to PRIVE2D%ADR(1)%P%R(I), despite it needing it to be 2D for now.
Everything works. Thanks again! |
|
The administrator has disabled public write access.
|
Assigning vertical head losses 7 years 7 months ago #28180
|
Hello!
Apologies for resurrecting this thread. Currently I'm increasing the number of horizontal layers in my mesh but when I specify 22 horizontal levels or more I get the "exceeding maximum iterations" error. Without the code it works fine. Decreasing the time step doesn't seem to do anything to resolve it. I can't figure out why changing to 22+ levels suddenly makes Telemac have a problem. Does anyone know why this might be? Best regards. |
|
The administrator has disabled public write access.
|
Assigning vertical head losses 7 years 7 months ago #28181
|
To add: previously I've had problems similar to this with mesh fidelity however this appears to be more linked to the fortran code. I also had a look at the mesh and adjusted the nodes in BlueKenue but no change either. The Selafin file is double precision also.
I've attached the relevant parts of my current code. It aims to ascribed head losses to the lower 4 layers of the mesh. Below is the example for a mesh with 34 horizontal layers. DO IPLAN = 2,5
TRANSF_PLANE%I(IPLAN)=3
ENDDO
ZPLANE%R(2)=-0.01D0
ZPLANE%R(3)=-0.01D0
ZPLANE%R(4)=-0.01D0
ZPLANE%R(5)=-0.01D0
DO IPLAN = 6,NPLAN
TRANSF_PLANE%I(IPLAN)=1
ENDDOS1U%TYPR='Q'
S1V%TYPR='Q'
S1W%TYPR='Q'
CALL OS( 'X=0 ',X=S1U)
CALL OS( 'X=0 ',X=S1V)
CALL OS( 'X=0 ',X=S1W)
IF(NONHYD) THEN
S0W%TYPR='0'
S1W%TYPR='Q'
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.6546
AREFV = 1.9889
AREFW = 13.8340
ELSE IF (IPLAN==(NPLAN-32))THEN
AREFU = 12.0348
AREFV = 15.0840
AREFW = 77.5191
ELSE IF (IPLAN==NPLAN-31) THEN
AREFU = 12.3361
AREFV = 14.8194
AREFW = 76.6161
ELSE IF (IPLAN==NPLAN-30) THEN
AREFU = 2.1330
AREFV = 2.4954
AREFW = 15.3330
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 DOBest regards. |
|
The administrator has disabled public write access.
|
Assigning vertical head losses 7 years 7 months ago #28183
|
Running debugger the error occurs in PREDIV, so it's something to do with how Telemac is trying to compute the free surface. This doesn't happen on all time-steps, but it always happens on the first time-step.
DANS NONHYDRO1 APPEL DE PREDIV GRACJG (BIEF) : EXCEEDING MAXIMUM ITERATIONS: 100 ABSOLUTE PRECISION: 0.1958074E-03 RETOUR DE PREDIV APPEL DE SOUKEP |
|
The administrator has disabled public write access.
|
|
Moderators: pham
