Friday, June 24, 2011

Programming the Initial Conditions in Telemac

The current project that I am working in involves looking at the vertical velocity profiles along a dune. By default Telemac3D creates equally spaced vertical layers from your model. This is not really desirable as  one generally wants a greater density of vertical layers near the bottom where the velocity profile changes rapidly and fewer near the free surface to save computational time. Telemac3D does allow you to specify the vertical profiles however you must do so programmatically - by creating a princi.f file and writing your own CONDMIN subroutine to initialize the initial conditions. After playing around with this I found a good example of the princi.f file with the examples and was able to modify it as required.

TELEMAC 3D uses a sigma transformation for the z or vertical access which means that at each point on the mesh z* = z/Zmax. This means that we specify our height from the bottom using values between 0 and 1.

!-----------------------------------------------------------------------
!
!     INITIALISATION OF ZSTAR, VERTICAL COORDINATE
!     IN SIGMA TRANSFORMATION
!
!     CASE WITHOUT REFERENCE PLANE
!     ----------------------------
!
!         ONE MUST HAVE :
!
!            * ZSTAR%R(1)     = 0.D0 ( BOTTOM PLANE )
!            * ZSTAR%R(NPLAN) = 1.D0 ( FREE SURFACE PLANE )
!
!         AND FOR ALL I BETWEEN 1 AND NPLAN-1
!
!            * ZSTAR%R(I) < ZSTAR%R(I+1)
!
!     CASE WITH REFERENCE PLANE
!     -------------------------
!
!         ONE MUST HAVE :
!
!            * ZSTAR%R(1)      = -1.D0 ( BOTTOM PLANE )
!            * ZSTAR%R(NPLINT) =  0.D0 ( REFERENCE PLANE )
!            * ZSTAR%R(NPLAN)  =  1.D0 ( FREE SURFACE PLANE )
!
!         AND FOR ALL I BETWEEN 1 AND NPLAN-1
!
!            * ZSTAR%R(I) < ZSTAR%R(I+1)
!
!     STANDARD IS: EVENLY SPACED PLANES
!
!     IF NO REFERENCE PLANE : NPLINT = 1
!
!     TRANSF IS KEY-WORD "MESH TRANSFORMATION"


 WRITE(lu,*) 'TRANSF IS KEY-WORD "MESH TRANSFORMATION"'
!
      IF(TRANSF.EQ.2) THEN
        WRITE(LU,*)
        IF(LNG.EQ.1) THEN
          WRITE(LU,*) 'AVEC TRANSFORMATION DU MAILLAGE = 2'
          WRITE(LU,*) 'DONNER DANS CONDIM LES VALEURS DE ZSTAR'
        ENDIF
        IF(LNG.EQ.2) THEN
          WRITE(LU,*) 'WITH MESH TRANSFORMATION = 2'
          WRITE(LU,*) 'GIVE THE VALUES OF ZSTAR IN CONDIM'
        ENDIF
        WRITE(LU,*)
!       CALL PLANTE(1)
!        STOP
!       EXAMPLE WITH 3 PLANES     
        ZSTAR%R(1)=0.D0
        ZSTAR%R(2)=0.01D0
        ZSTAR%R(3)=0.03D0 
        ZSTAR%R(4)=0.05D0
        ZSTAR%R(5)=0.07D0
        ZSTAR%R(6)=0.09D0
        ZSTAR%R(7)=0.11D0
        ZSTAR%R(8)=0.13D0 
        ZSTAR%R(9)=0.15D0
        ZSTAR%R(10)=0.25D0
        ZSTAR%R(11)=0.368D0
        ZSTAR%R(12)=0.5D0
        ZSTAR%R(13)=0.75D0
        ZSTAR%R(14)=1.D0
      ELSEIF(TRANSF.EQ.3) THEN
        IF(NPLINT.GE.2) THEN
          DO IPLAN = 1,NPLINT-1
            ZSTAR%R(IPLAN) = DBLE(IPLAN-NPLINT)/DBLE(NPLINT-1)
          ENDDO
        ENDIF
        DO IPLAN = NPLINT,NPLAN
          ZSTAR%R(IPLAN) = DBLE(IPLAN-NPLINT)/DBLE(NPLAN-NPLINT)
        ENDDO
      ELSEIF(TRANSF.EQ.4) THEN
        WRITE(LU,*)
        IF(LNG.EQ.1) THEN
          WRITE(LU,*) 'AVEC TRANSFORMATION DU MAILLAGE = 4'
          WRITE(LU,*) 'DONNER DANS CONDIM LES VALEURS DE ZSTAR'
          WRITE(LU,*) 'COTES REELLES SOOUHAITEES DES PLANS'
        ENDIF
        IF(LNG.EQ.2) THEN
          WRITE(LU,*) 'WITH MESH TRANSFORMATION = 4'
          WRITE(LU,*) 'GIVE THE VALUES OF ZSTAR IN CONDIM'
          WRITE(LU,*) 'REAL ELEVATION OF PLANES'
        ENDIF
        WRITE(LU,*)
        CALL PLANTE(1)
        STOP
!        EXAMPLE WITH 4 PLANES        
!        ZSTAR%R(1)=-10.D0
!        ZSTAR%R(2)=-9.D0
!        ZSTAR%R(3)=-8.D0
!        ZSTAR%R(4)=-7.D0
      ELSE
!       SIGMA TRANSFORMATION AND OTHER CASES
        DO IPLAN = 1,NPLAN
          ZSTAR%R(IPLAN) = DBLE(IPLAN-1)/DBLE(NPLAN-1)
        ENDDO
      ENDIF

No comments:

Post a Comment