bnd_spectra.f

Go to the documentation of this file.

!                   ******************
                    MODULE bnd_spectra
!                   ******************
!
!
!***********************************************************************
!     TOMAWAC   V7P3                                   23/02/2017
!***********************************************************************
!
!brief    MODULE TO IMPOSE SPECTRA ON OPEN BOUNDARIES FROM AN
!         EXTERNAL MESH FILE
!
!     history  A. JOLY (EDF - LNHE)
!     +        23/02/2017
!     +        V7P3
!     +   CREATED
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!
      IMPLICIT NONE
!
      PRIVATE
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      PUBLIC :: impose_bnd_spectra
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      ! TIME VARIABLES
      DOUBLE PRECISION :: tv1,tv2
      INTEGER :: record1,record2
      DOUBLE PRECISION :: time1,time2
      !ERROR INTEGER WHEN READING THE MESH
      INTEGER :: ierr
      ! BOUNDARY SPECTRA VARIABLES
      CHARACTER(LEN=80) :: bc_title ! TITLE
      INTEGER :: bc_nvar   ! NUMBER OF VARIABLES
      INTEGER :: bc_npoin  ! NUMBER OF MESH NODES
      INTEGER :: bc_typ_elem  ! TYPE OF ELEMENT
      INTEGER :: bc_nelem  ! NUMBER OF ELEMENTS
      INTEGER :: bc_ndp    ! NUMBER OF ELEMENT FACES
      INTEGER :: bc_nptfr  ! NUMBER OF BOUNDARY NODES
      INTEGER :: bc_nptir  ! NUMBER OF INTERFACES
      INTEGER :: bc_ndim  ! DIMENSION OF THE MESH
      CHARACTER(LEN=16), ALLOCATABLE :: bc_varlist(:)
      CHARACTER(LEN=16), ALLOCATABLE :: bc_unitlist(:)
      INTEGER :: bc_nt  ! NUMBER TIME STEPS
      DOUBLE PRECISION, ALLOCATABLE :: bc_x(:),bc_y(:) !X AND Y COORDINATES ON SPE MESH
      DOUBLE PRECISION, ALLOCATABLE :: bc_freq(:),bc_teta(:) !F AND THETA COORDINATES OF SPECTRUM
      DOUBLE PRECISION, ALLOCATABLE :: bc_val(:)  ! TEMPORY DATA FOR READING
      DOUBLE PRECISION, ALLOCATABLE :: bc_all1(:,:), bc_all2(:,:) !VALUES A TIME TV1 AND TV2
      !FOR POINT OF THE SPECTRA FILE
      INTEGER, ALLOCATABLE :: bc_corr(:) ! INDEX OF THE CLOSEST SPECTRUM TO A LIQUID BOUNDARY NODE
      ! OPTIONAL ARGUMENTS
      INTEGER nnelebd, typ_bnd
!

      INTEGER :: nptfr_loc
      INTEGER, ALLOCATABLE :: ind_ptfr(:) ! INDEX INTO VARLIST AND BOUNDARY NUMBER
!
      SAVE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
!***********************************************************************
      CONTAINS
!***********************************************************************

!                   *****************************
                    SUBROUTINE impose_bnd_spectra
!                   *****************************
!
     &     (imp_file,lt,at,fbor,nptfr,ndire,nf)
!
!***********************************************************************
!     TOMAWAC   V7P3                                   23/02/2017
!***********************************************************************
!
!     brief    READS SPECTRA SAVED IN A MESH FILE AND IMPOSES THEM ON
!     +        OPEN BOUNDARY POINTS
!
!     history  A. JOLY (EDF - LNHE)
!     +        23/02/2017
!     +        V7P3
!     +   CREATED
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| IMP_FILE       |-->| MESH FILE WITH THE IMPOSED SPECTRA
!| LT             |-->| NUMBER OF THE TIME STEP CURRENTLY SOLVED
!| AT             |-->| COMPUTATION TIME
!| FBOR           |<->| SPECTRAL VARIANCE DENSITY AT THE BOUNDARIES
!| NPTFR          |-->| NUMBER OF BOUNDARY POINTS
!| NDIRE          |-->| NUMBER OF DIRECTIONS
!| NF             |-->| NUMBER OF FREQUENCIES
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE interface_hermes
      USE declarations_special
      USE declarations_telemac
      USE bief_def, ONLY : bief_file
      USE declarations_tomawac, ONLY : npspe,xspe,yspe,unitspe,phasspe,
     &                                 x,y,nbor,lifbor,f1,raisf
!
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      TYPE(bief_file), INTENT(IN)    :: IMP_FILE
      INTEGER, INTENT(IN)            :: LT
      DOUBLE PRECISION, INTENT(IN)   :: AT
      INTEGER, INTENT(IN)            :: NPTFR,NDIRE,NF
      DOUBLE PRECISION, INTENT(INOUT):: FBOR(nptfr,ndire,nf)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER IPTFR, ISPE, IPOIN, IDIR, IFF,ITMP
      INTEGER X_ORIG, Y_ORIG
      DOUBLE PRECISION DIST_SPE,DIST
      DOUBLE PRECISION COEFT
      DOUBLE PRECISION FCL1,FCL2
      INTEGER BC_NF,BC_NDIRE,ITRAV
      DOUBLE PRECISION BC_F1,BC_F2,BC_RAISF
      DOUBLE PRECISION F_TEMP
      DOUBLE PRECISION, PARAMETER :: EPS = 1.d-6
!
!-----------------------------------------------------------------------
! DURING THE FIRST TIME STEP
!-----------------------------------------------------------------------
      IF(lt.EQ.0)THEN
! 1/ FIND THE CLOSEST SPECTRUM FOR EACH LIQUID BOUNDARY NODE
        IF(ALLOCATED(bc_corr))DEALLOCATE(bc_corr)
        ALLOCATE(bc_corr(nptfr))
        DO iptfr=1,nptfr
          IF(lifbor(iptfr).EQ.kent)THEN
            dist_spe=1.d99
            DO ispe=1,npspe
              dist=sqrt((xspe(ispe)-x(nbor(iptfr)))**2+
     &                  (yspe(ispe)-y(nbor(iptfr)))**2)
              IF(dist.LE.dist_spe)THEN
                bc_corr(iptfr)=ispe
                dist_spe=dist
              ENDIF
            ENDDO
          ELSE
            bc_corr(iptfr)=0
          ENDIF
        END DO
! 2/ READ THE MESH THE SPECTRA TO BE IMPOSED ON THE BOUNDARY
        ! MESH INFO
        CALL read_mesh_info(imp_file%FMT,imp_file%LU,
     &          bc_title,bc_nvar,
     &          bc_npoin,bc_typ_elem,bc_nelem,bc_nptfr,bc_nptir,
     &          bc_ndp,itrav,x_orig,y_orig,typ_bnd,nnelebd)
        CALL get_mesh_dimension(imp_file%FMT,imp_file%LU,
     &          bc_ndim,ierr)
        ! READ THE NAME OF THE VARIABLES
        IF(ALLOCATED(bc_varlist))DEALLOCATE(bc_varlist)
        IF(ALLOCATED(bc_unitlist))DEALLOCATE(bc_unitlist)
        ALLOCATE(bc_varlist(bc_nvar))
        ALLOCATE(bc_unitlist(bc_nvar))
        CALL get_data_var_list(imp_file%FMT,imp_file%LU,
     &          bc_nvar,bc_varlist,bc_unitlist,ierr)
! 3/ READ THE NUMBER OF TIME STEPS AND THE FIRST TIME
        ! GET NUMBER OF TIME STEPS
        CALL get_data_ntimestep(imp_file%FMT,imp_file%LU,
     &          bc_nt,ierr)
        ! READ THE FIRST TIME
        record1 = 0
        record2 = 1
        CALL get_data_time(imp_file%FMT,imp_file%LU,
     &          record1,time1,ierr)
        tv1=(time1-phasspe)*unitspe

        IF(tv1.GT.at) THEN
          WRITE(lu,*)'*********************************************'
          WRITE(lu,*)'THE FIRST RECORDING OF THE FILE'
          WRITE(lu,*)'  ',tv1,' IS OLDER THAN THE BEGINNING'
          WRITE(lu,*)'  OF THE COMPUTATION',at
          WRITE(lu,*)'*********************************************'
          CALL plante(1)
          stop
        ENDIF

        DO
        CALL get_data_time(imp_file%FMT,imp_file%LU,
     &          record2,time2,ierr)
        tv2=(time2-phasspe)*unitspe
        IF(tv2.LT.at) THEN
          record1 = record2
          record2 = record2 + 1
          tv1 = tv2
          IF (record2.GT.bc_nt) THEN
            WRITE(lu,*)'*****************************************'
            IF(lng.EQ.1) THEN
              WRITE(lu,*)'LA FIN DU FICHIER DE CONDITION LIMITE'
              WRITE(lu,*)'EST ATTEINTE AU TEMPS : ',at
            ELSE
              WRITE(lu,*)'THE END OF THE BOUNDARY CONDITION FILE'
              WRITE(lu,*)'IS REACHED AT TIME: ',at
            ENDIF
            WRITE(lu,*)'*****************************************'
            CALL plante(1)
            stop
          ENDIF
        ELSE
          EXIT
        ENDIF
      ENDDO


! 4/ READ X AND Y COORDINATES
        IF(ALLOCATED(bc_x))DEALLOCATE(bc_x)
        IF(ALLOCATED(bc_y))DEALLOCATE(bc_y)
        IF(ALLOCATED(bc_freq))DEALLOCATE(bc_freq)
        IF(ALLOCATED(bc_teta))DEALLOCATE(bc_teta)
        ALLOCATE(bc_x(bc_npoin))
        ALLOCATE(bc_y(bc_npoin))
        ALLOCATE(bc_freq(bc_npoin))
        ALLOCATE(bc_teta(bc_npoin))
!
        CALL get_mesh_coord(imp_file%FMT,imp_file%LU,
     &            1,bc_ndim,bc_npoin,bc_x,ierr)
        CALL get_mesh_coord(imp_file%FMT,imp_file%LU,
     &          2,bc_ndim,bc_npoin,bc_y,ierr)
! 5/ CHECK THAT THE MESH READ HAS THE SAME THE FREQUENCIES AND
!    DIRECTIONS AS IN THE SIMULATION
        bc_nf=0
        bc_ndire=0
        bc_f1=10.d10
        bc_f2=10.d10
        bc_raisf=0.d0
        DO ipoin=1,bc_npoin
          IF((abs(bc_x(ipoin)).LE.eps).AND.(bc_y(ipoin).GE.0.0))THEN
            bc_nf=bc_nf+1
            f_temp=bc_y(ipoin)
            IF(f_temp.LT.bc_f1)THEN
              bc_f2=bc_f1
              bc_f1=f_temp
            ELSEIF(f_temp.LT.bc_f2)THEN
              bc_f2=f_temp
            ENDIF
          ENDIF
        END DO
        bc_raisf=bc_f2/bc_f1
        bc_ndire=bc_npoin/bc_nf
        IF(abs(bc_f1-f1).GT.1.d-6) THEN
          WRITE(lu,*)'*********************************************'
          WRITE(lu,*)'THE MINIMAL FREQUENCY OF THE SPECTRA READ'
          WRITE(lu,*)'IS NOT EQUAL TO THE FREQUENCY OF THE'
          WRITE(lu,*)'SIMULATION.'
          WRITE(lu,*)'- FREQUENCY OF THE SIMULATION:',f1
          WRITE(lu,*)'- FREQUENCY READ:',bc_f1
          WRITE(lu,*)'*********************************************'
          CALL plante(1)
          stop
        ENDIF
        IF(abs(bc_raisf-raisf).GT.1.d-6) THEN
          WRITE(lu,*)'*********************************************'
          WRITE(lu,*)'THE FREQUENTIAL RATIO OF THE SPECTRA READ'
          WRITE(lu,*)'IS NOT EQUAL TO THE RATIO OF THE SIMULATION.'
          WRITE(lu,*)'- FREQUENTIAL RATIO OF THE SIMULATION:',raisf
          WRITE(lu,*)'- FREQUENTIAL RATIO READ:',bc_raisf
          WRITE(lu,*)'*********************************************'
          CALL plante(1)
          stop
        ENDIF
        IF(bc_nf.NE.nf) THEN
          WRITE(lu,*)'*********************************************'
          WRITE(lu,*)'THE NUMBER OF FREQUENCY OF THE SPECTRA READ'
          WRITE(lu,*)'IS NOT EQUAL TO THE NUMBER OF FREQUENCY OF THE'
          WRITE(lu,*)'SIMULATION.'
          WRITE(lu,*)'- NUMBER OF FREQUENCY OF THE SIMULATION:',nf
          WRITE(lu,*)'- NUMBER OF FREQUENCY READ:',bc_nf
          WRITE(lu,*)'*********************************************'
          CALL plante(1)
          stop
        ENDIF
        IF(bc_ndire.NE.ndire) THEN
          WRITE(lu,*)'*********************************************'
          WRITE(lu,*)'THE NUMBER OF DIRECTIONS OF THE SPECTRA READ'
          WRITE(lu,*)'IS NOT EQUAL TO THE NUMBER OF DIRECTIONS OF THE'
          WRITE(lu,*)'SIMULATION.'
          WRITE(lu,*)'- NUMBER OF DIRECTIONS OF THE SIMULATION:',ndire
          WRITE(lu,*)'- NUMBER OF DIRECTIONS READ:',bc_ndire
          WRITE(lu,*)'*********************************************'
          CALL plante(1)
          stop
        ENDIF
! 6/ CALCULATE FREQ AND DIR
        DO ipoin=1,bc_npoin
          bc_teta(ipoin)=atan2(bc_x(ipoin),bc_y(ipoin))+atan(1.d0)*8.d0
          bc_teta(ipoin)=mod(bc_teta(ipoin),atan(1.d0)*8.d0)
          IF(sin(bc_teta(ipoin)).NE.0.d0)THEN
            bc_freq(ipoin)=bc_x(ipoin)/sin(bc_teta(ipoin))
          ELSE
            bc_freq(ipoin)=bc_y(ipoin)/cos(bc_teta(ipoin))
          ENDIF
        ENDDO
! 6.5/ CHECK NUMBER OF POINTS ON PROCESS
        nptfr_loc  = 0
        ALLOCATE(ind_ptfr(nptfr))
        DO iptfr=1,nptfr
          IF(lifbor(iptfr).EQ.kent)THEN
            nptfr_loc = nptfr_loc+1
            ind_ptfr(nptfr_loc) = iptfr
          ENDIF
        ENDDO


! 7/ ALLOCATE THE TEMPORARY SPECTRUM AT THE TWO TIME STEPS
        IF(ALLOCATED(bc_val))DEALLOCATE(bc_val)
        ALLOCATE(bc_val(bc_npoin))

        IF(ALLOCATED(bc_all1))DEALLOCATE(bc_all1)
        IF(ALLOCATED(bc_all2))DEALLOCATE(bc_all2)
        ALLOCATE(bc_all1(bc_npoin,nptfr_loc))
        ALLOCATE(bc_all2(bc_npoin,nptfr_loc))

! 8/ READ FIRST TIME
        DO iptfr=1,nptfr_loc
          itmp = ind_ptfr(iptfr)
          ! READ DATA BEFORE
          CALL get_data_value(imp_file%FMT,imp_file%LU,
     &       record1,bc_varlist(bc_corr(itmp)),bc_val,bc_npoin,ierr)
          DO ipoin = 1,bc_npoin
            bc_all1(ipoin,iptfr)= bc_val(ipoin)
          ENDDO
          ! READ DATA AFTER
          CALL get_data_value(imp_file%FMT,imp_file%LU,
     &       record2,bc_varlist(bc_corr(itmp)),bc_val,bc_npoin,ierr)
          DO ipoin = 1,bc_npoin
            bc_all2(ipoin,iptfr)= bc_val(ipoin)
          ENDDO
        ENDDO


!-----------------------------------------------------------------------
      ENDIF
!-----------------------------------------------------------------------
! FOR ALL TIME STEPS
!-----------------------------------------------------------------------
! 8/ READ THE SECOND TIME STEP

      ! RETURN IF NO BOUNDARY DATA
      IF (nptfr_loc.EQ.0) THEN
            RETURN
      ENDIF

! 8.5/ UPDATE DATA IF NEEDED
      IF (at.GT.tv2) THEN
        DO
          record2 = record2+1
          tv1 = tv2
          CALL get_data_time(imp_file%FMT,imp_file%LU,
     &            record2,time2,ierr)
          tv2=(time2-phasspe)*unitspe
          IF(tv2.LT.at) THEN
            IF (record2.GT.bc_nt) THEN
              WRITE(lu,*)'*****************************************'
              IF(lng.EQ.1) THEN
                WRITE(lu,*)'LA FIN DU FICHIER DE CONDITION LIMITE'
                WRITE(lu,*)'EST ATTEINTE AU TEMPS : ',at
              ELSE
                WRITE(lu,*)'THE END OF THE BOUNDARY CONDITION FILE'
                WRITE(lu,*)'IS REACHED AT TIME: ',at
              ENDIF
              WRITE(lu,*)'*****************************************'
              CALL plante(1)
              stop
            ENDIF
          ELSE
            EXIT
          ENDIF
        ENDDO

        IF (record1.LT.record2-1) THEN
            ! IN CASE TIME STEP IN SPECTRAL FILE IS SMALLER THAN MODEL TIME STEP
            ! ALSO UPDATE RECORD 1
            record1 = record2-1
            DO iptfr=1,nptfr_loc
              itmp = ind_ptfr(iptfr)
              CALL get_data_value(imp_file%FMT,imp_file%LU,
     &           record1,bc_varlist(bc_corr(itmp)),bc_val,
     &           bc_npoin,ierr)
              DO ipoin = 1,bc_npoin
                bc_all1(ipoin,iptfr)= bc_val(ipoin)
              ENDDO
            ENDDO
        ELSE
          ! COPY DATA
          record1 = record2
          DO ipoin = 1,bc_npoin
            bc_all2(ipoin,iptfr)= bc_val(ipoin)
          ENDDO
        ENDIF
        ! UPDATE NEW TIME
        DO iptfr=1,nptfr_loc
          itmp = ind_ptfr(iptfr)
          CALL get_data_value(imp_file%FMT,imp_file%LU,
     &       record2,bc_varlist(bc_corr(itmp)),bc_val,bc_npoin,ierr)
          DO ipoin = 1,bc_npoin
            bc_all2(ipoin,iptfr)= bc_val(ipoin)
          ENDDO
        ENDDO
      ENDIF
!
      coeft=(at-tv1)/(tv2-tv1)
! 9/ IMPOSE THE SPECTRA
        !TODO; USE BETTER INTERPOLATION
      DO iptfr=1,nptfr_loc
          itmp = ind_ptfr(iptfr)
          DO idir=1,ndire
            DO iff=1,nf
              ipoin=(idir+ndire*(iff-1))
              fcl1=bc_all1(ipoin,iptfr)
              fcl2=bc_all2(ipoin,iptfr)
              fbor(itmp,idir,iff)=fcl1+(fcl2-fcl1)*coeft
            ENDDO
          ENDDO
      ENDDO

!
      RETURN
      END SUBROUTINE impose_bnd_spectra
!
      END MODULE bnd_spectra