gretel_autop.f

Go to the documentation of this file.

!                   ***********************
                    SUBROUTINE gretel_autop
!                   ***********************
     &(geo,geoformat,bnd,res,resformat,nproc,nplan_res,method)
!
!
!***********************************************************************
! PARALLEL   V6P2                                   21/08/2010
!***********************************************************************
!
!brief    MERGES THE RESULTS OF A PARALLEL COMPUTATION
!+                TO WRITE A SINGLE FILE IN A GIVEN FORMAT.
!
!
!history  Y. Audouin
!+        02/09/2014
!+
!+   Creation of the file
!+   This version of gretel now handles multiple format and has a better
!+   memory organisation
!+   It also uses the hermes module for I/O
!
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE interface_hermes
      USE declarations_special
      USE bief, ONLY: read_mesh_info
!
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
!
      CHARACTER(LEN=PATH_LEN), INTENT(IN) :: GEO
      CHARACTER(LEN=PATH_LEN), INTENT(IN) :: BND
      CHARACTER(LEN=PATH_LEN), INTENT(IN) :: RES
      CHARACTER(LEN=8),   INTENT(INOUT) :: GEOFORMAT,RESFORMAT
      INTEGER,            INTENT(IN) :: NPROC
      INTEGER,            INTENT(INOUT) :: NPLAN_RES
      INTEGER,            INTENT(IN) :: METHOD
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
!
      INTEGER IPID
      INTEGER I,J,IELEM
      INTEGER NPLAN_GEO,NELEM_GEO,NDP,NELEBD, NPTFR, NPTIR
      INTEGER :: NDIM, X_ORIG, Y_ORIG
!
      INTEGER, DIMENSION(:)  , ALLOCATABLE :: IPOBO_GEO,IPOBO3D
      INTEGER, DIMENSION(:), ALLOCATABLE :: KNOLG, TMP2
      INTEGER, DIMENSION(:), ALLOCATABLE :: IKLE_GEO,IKLE3D
      INTEGER, DIMENSION(:), ALLOCATABLE :: IKLE_BND
!
      DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: TMP, X, Y
!
      DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: X3D, Y3D, Z3D
!
      CHARACTER(LEN=300) :: RESPAR
!
      CHARACTER(LEN=80) TITSEL
      CHARACTER(LEN=32),ALLOCATABLE :: TEXTELU(:)
      CHARACTER(LEN=16),ALLOCATABLE :: VAR_NAME(:), VAR_UNIT(:)
      CHARACTER(LEN=16) :: VARNAME
      CHARACTER(LEN=11) EXTENS
      EXTERNAL  extens
      INTRINSIC real
!
      INTEGER IERR, NRES, NRESPAR, NGEO
      INTEGER TYP_ELEM, TYP_BND_ELEM
      INTEGER DATE(3), TIME(3)
      INTEGER NTIMESTEP_RES
      INTEGER NPOIN_GEO, NPOIN_RES, NPOIN_PAR
      INTEGER NVAR_RES,NVAR_GEO
      INTEGER NPOIN3D, NELEM3D
      INTEGER DATE_TMP(6)
!
!-------------------------------------------------------------------------
!
!
!|==================================================================|
!|                                                                  |
!| START: MERGES FILES RESULTING FROM THE DOMAIN DECOMPOSITION      |
!|                                                                  |
!|==================================================================|
!
! READS FILE NAMES AND THE NUMBER OF PROCESSORS / PARTITIONS
!
!
!
!     Header information
!
      respar = trim(res) // extens(nproc-1,0)
      CALL open_mesh(resformat,respar,nrespar,'READ     ',ierr)
      CALL check_call(ierr,"GRETEL:OPEN_MESH:RESPAR")
!     Opening result file after so that it will be the same endianess as
!     the partitionned result files
      CALL open_mesh(resformat,res,nres,'WRITE    ',ierr)
      CALL check_call(ierr,"GRETEL:OPEN_MESH:RES")
!
      ! GET THE MESH TITLE
      CALL get_mesh_title(resformat,nrespar,titsel,ierr)
      CALL check_call(ierr,"GRETEL:GET_MESH_TITLE")

!     Get the number of variable
      CALL get_data_nvar(resformat,nrespar,nvar_res,ierr)
      CALL check_call(ierr,"GRETEL:GET_DATA_NVAR")
!
      ALLOCATE(textelu(nvar_res),stat=ierr)
      CALL check_allocate(ierr,'GRETEL:TEXTELU')
      ALLOCATE(var_name(nvar_res),stat=ierr)
      CALL check_allocate(ierr,'GRETEL:VAR_NAME')
      ALLOCATE(var_unit(nvar_res),stat=ierr)
      CALL check_allocate(ierr,'GRETEL:VAR_UNIT')
!
      CALL get_data_var_list(resformat,nrespar,nvar_res,var_name,
     &                       var_unit,ierr)
      CALL check_call(ierr,"GRETEL:GET_DATA_VAR_LIST")
!
      WRITE(lu,*) 'TITLE=',titsel
      WRITE(lu,*) 'NBVAR=',nvar_res
      DO i=1,nvar_res
        textelu(i)(1:16) = var_name(i)
        textelu(i)(17:32) = var_unit(i)
        WRITE(lu,*) 'VARIABLE ',i,' : ',textelu(i)
      ENDDO ! I
!
      DEALLOCATE(var_name)
      DEALLOCATE(var_unit)
!
!     WE NEED TO GET THE NUMBER OF PLANES IN THE PARTITIONNED FILE
!     TO KNOW IF WE NEED TO TRANFORM THE GEOMETRY IN 3D GEOMETRY
      CALL get_mesh_nplan(resformat,nrespar,nplan_res,ierr)
      CALL check_call(ierr,"GRETEL:MESH_NPLAN")
!
!     Get the number of timestep
      CALL get_data_ntimestep(resformat,nrespar,ntimestep_res,ierr)
      CALL check_call(ierr,"GRETEL:GET_DATA_NTIMESTEP")

      CALL close_mesh(resformat,nrespar,ierr)
      CALL check_call(ierr,"GRETEL:CLOSE_MESH:RESPAR")

      CALL set_header(resformat,nres,titsel,nvar_res,textelu,ierr)
!
!     Geometry information
!
      CALL open_mesh(geoformat,geo,ngeo,'READ     ',ierr)
      CALL check_call(ierr,"GRETEL:OPEN_MESH:GEO")
!

!     Boundary file
      CALL open_bnd(geoformat,bnd,ngeo,'READ     ',ierr)
      CALL check_call(ierr,"GRETEL:OPEN_BND:GEO")

      CALL read_mesh_info(geoformat,ngeo,titsel,nvar_geo,npoin_geo,
     &                    typ_elem,nelem_geo,nptfr,nptir,ndp,nplan_geo,
     &                    x_orig,y_orig,typ_bnd_elem,nelebd)
!
      WRITE(lu,*) 'GEO MESH INFORMATIONS:'
      WRITE(lu,*) 'NELEM=',nelem_geo
      WRITE(lu,*) 'NPOIN=',npoin_geo
      WRITE(lu,*) 'NDP=',ndp
      WRITE(lu,*) 'TYP_ELEM=',typ_elem
!
      ALLOCATE(ikle_geo(nelem_geo*ndp),stat=ierr)
      CALL check_allocate(ierr,'GRETEL:IKLE_GEO')
      ALLOCATE(ipobo_geo(npoin_geo),stat=ierr)
      CALL check_allocate(ierr,'GRETEL:IPOBO')
!
      CALL read_mesh_conn(geoformat,ngeo,npoin_geo,typ_elem,nelem_geo,
     &                    ndp,typ_bnd_elem,nelebd,ikle_geo, ipobo_geo)

      IF(nplan_res.LE.1) THEN
        ndim = 2
      ELSE
        ndim = 3
      ENDIF
      ! GET 2D COORDINATES
      ALLOCATE(x(npoin_geo),stat=ierr)
      CALL check_allocate(ierr,'GRETEL:X')
      ALLOCATE(y(npoin_geo),stat=ierr)
      CALL check_allocate(ierr,'GRETEL:Y')
      ALLOCATE(z3d(npoin_geo*nplan_res),stat=ierr)
      CALL check_allocate(ierr,'GRETEL:Z3D')
      ! GET MESH COORDINATES FROM THE GEO MESH
      CALL get_mesh_coord(geoformat,ngeo,1,2,npoin_geo,x,ierr)
      CALL check_call(ierr,'GRETEL:GET_MESH_COORD:X:GEO')
      CALL get_mesh_coord(geoformat,ngeo,2,2,npoin_geo,y,ierr)
      CALL check_call(ierr,'GRETEL:GET_MESH_COORD:Y:GEO')
      !
      ! Update coordiantes with coordinates from the partitionned files
      ! as they could have been modified by corrxy
      DO ipid = 0, nproc-1
!
        respar = trim(res) // extens(nproc-1,ipid)
        CALL open_mesh(resformat,respar,nrespar,'READ     ',ierr)
        CALL check_call(ierr,"GRETEL:OPEN_MESH:RESPAR2")
!
        CALL get_mesh_npoin(resformat,nrespar,typ_elem,npoin_par,ierr)
        CALL check_call(ierr,"GRETEL:GET_MESH_NPOIN:RESPAR")
!
        ALLOCATE(knolg(npoin_par),stat=ierr)
        CALL check_allocate(ierr,'GRETEL:KNOLG')
        ALLOCATE(tmp(npoin_par),stat=ierr)
        CALL check_allocate(ierr,'GRETEL:TMP')
!
        CALL get_mesh_l2g_numbering(resformat,nrespar,knolg,
     &                              npoin_par,ierr)
        CALL check_call(ierr,'GRETEL:GET_MESH_L2G_NUMBERING:RESPAR')
!
        CALL get_mesh_coord(resformat,nrespar,1,ndim,npoin_par,tmp,ierr)
        CALL check_call(ierr,'GRETEL:GET_MESH_COORD:X:RESPAR')
        DO i=1,npoin_par/(max(nplan_res,1))
          x(knolg(i)) = tmp(i)
        ENDDO
        CALL get_mesh_coord(resformat,nrespar,2,ndim,npoin_par,tmp,ierr)
        CALL check_call(ierr,'GRETEL:GET_MESH_COORD:Y:RESPAR')
        DO i=1,npoin_par/(max(nplan_res,1))
          y(knolg(i)) = tmp(i)
        ENDDO
        ! If 3d getting Z from result file (First variable at time=0.0)
        IF(ndim.EQ.3) THEN
          varname = textelu(1)(1:16)
          CALL get_data_value(resformat,nrespar,0,varname,tmp,
     &                        npoin_par,ierr)
          DO i=1,npoin_par
            z3d(knolg(i)) = tmp(i)
          ENDDO
        ENDIF
        ! Getting the date from result file
        CALL get_mesh_date(resformat,nrespar,date_tmp,ierr)
        CALL check_call(ierr,'GRETEL:GET_MESH_DATE;RESPAR')
        DO i=1,3
          date(i) = date_tmp(i)
          time(i) = date_tmp(i+3)
        ENDDO
        !
        CALL close_mesh(resformat,nrespar,ierr)
        CALL check_call(ierr,'GRETEL:CLOSEMESH:RESPAR')
        DEALLOCATE(tmp)
        DEALLOCATE(knolg)
      ENDDO ! IPID
      !
      ! IF WE HAVE A 3D RESULT WE NEED TO TRANSFORM THE MESH IN 3D
      ! WRITES THE MESH INFORMATION TO THE MERGED FILE
      WRITE(lu,*) 'WRITING MESH'
      IF(ndim.EQ.2) THEN
        ! 2D
        ALLOCATE(tmp2(nelem_geo*ndp),stat=ierr)
        CALL check_allocate(ierr,'GRETEL:TMP0')
        DO i = 1,ndp
          DO ielem = 1,nelem_geo
            tmp2((i-1)*nelem_geo + ielem) = ikle_geo((ielem-1)*ndp+i)
          ENDDO
        ENDDO

        CALL set_mesh(resformat,nres,2,typ_elem,ndp,nptfr,nptir,
     &                nelem_geo,npoin_geo,tmp2,ipobo_geo,ipobo_geo,x,y,
     &                nplan_res,date,time,x_orig,y_orig,ierr,
     &                in_place=.true.)
        CALL check_call(ierr,'GRETEL:SET_MESH:RES')
        DEALLOCATE(ikle_geo)
        DEALLOCATE(tmp2)
        DEALLOCATE(ipobo_geo)
        DEALLOCATE(x)
        DEALLOCATE(y)
        npoin_res = npoin_geo
      ELSE
        ! 3D
        ndp = 6
        typ_elem = prism_elt_type
        npoin3d = npoin_geo * nplan_res
        nelem3d = nelem_geo * (nplan_res-1)
        nptfr = 0
        nptir = 0
        npoin_res = npoin3d
        ALLOCATE(ikle3d(nelem3d*ndp),stat=ierr)
        CALL check_allocate(ierr, 'GRETEL:IKLE3D')
        ALLOCATE(ipobo3d(npoin3d),stat=ierr)
        CALL check_allocate(ierr, 'GRETEL:IPOBO3D')
!
        ! BUILDING 3D IKLE
        DO i=1,nplan_res-1
          DO j=1,nelem_geo
            ikle3d(j+(i-1)*nelem_geo + (0*nelem3d)) =
     &                    ikle_geo((j-1)*3+1) + (i-1)*npoin_geo
            ikle3d(j+(i-1)*nelem_geo + (1*nelem3d)) =
     &                    ikle_geo((j-1)*3+2) + (i-1)*npoin_geo
            ikle3d(j+(i-1)*nelem_geo + (2*nelem3d)) =
     &                    ikle_geo((j-1)*3+3) + (i-1)*npoin_geo
            ikle3d(j+(i-1)*nelem_geo + (3*nelem3d)) =
     &                    ikle_geo((j-1)*3+1) + (i)*npoin_geo
            ikle3d(j+(i-1)*nelem_geo + (4*nelem3d)) =
     &                    ikle_geo((j-1)*3+2) + (i)*npoin_geo
            ikle3d(j+(i-1)*nelem_geo + (5*nelem3d)) =
     &                    ikle_geo((j-1)*3+3) + (i)*npoin_geo
          ENDDO
        ENDDO
        ! DEFAULT IPOBO
        DO i=1,npoin3d
          ipobo3d(i) = 0
        ENDDO
        ! WE DONT NEED THEM ANYMORE
        DEALLOCATE(ikle_geo)
        DEALLOCATE(ipobo_geo)
        ! BUILD THE COORDINATES FROM THE 2D COORDINATES
        ALLOCATE(x3d(npoin_geo*nplan_res),stat=ierr)
        CALL check_allocate(ierr,'GRETEL:X3D')
        ALLOCATE(y3d(npoin_geo*nplan_res),stat=ierr)
        CALL check_allocate(ierr,'GRETEL:Y3D')
        !
        DO i=1,npoin_geo
          DO j=1,nplan_res
            x3d(i+(j-1)*npoin_geo) = x(i)
            y3d(i+(j-1)*npoin_geo) = y(i)
          ENDDO
        ENDDO
        DEALLOCATE(x)
        DEALLOCATE(y)
        date = (/0,0,0/)
        time = (/0,0,0/)
        CALL set_mesh(resformat,nres,3,typ_elem,ndp,nptfr,nptir,nelem3d,
     &                npoin3d,ikle3d,ipobo3d,ipobo3d,x3d,y3d,
     &                nplan_res,date,time,x_orig,y_orig,ierr,
     &                z=z3d,in_place=.true.)
        CALL check_call(ierr,'GRETEL:SET_MESH:RES')
        DEALLOCATE(ikle3d)
        DEALLOCATE(ipobo3d)
        DEALLOCATE(x3d)
        DEALLOCATE(y3d)
        DEALLOCATE(z3d)


      ENDIF
      ! Transfering boundary information
      ! Getting boundary connectivity
      IF(geoformat.EQ.resformat) THEN
        ALLOCATE(ikle_bnd(nelebd*2), stat=ierr)
        CALL get_bnd_connectivity(geoformat, ngeo, typ_bnd_elem, nelebd,
     &                            2, ikle_bnd, ierr)
        CALL check_call(ierr,'GRETEL:GET_BND_CONNECTIVITY:GEO')

        CALL transfer_group_info(geoformat, ngeo, nres, typ_elem,
     &   typ_bnd_elem,ikle_bnd,nelebd,2,.true.,.true.,ierr)
        CALL check_call(ierr,'GRETEL:TRANSFER_GROUP_INFO:GEO')
        DEALLOCATE(ikle_bnd)
      ENDIF
!
!     Read results informations from partitioned files
!
      CALL merge_data
     &(npoin_res, nvar_res, ntimestep_res, nproc, resformat, nres,
     & typ_elem, textelu, res, ndim, nplan_res, npoin_geo, method)

      ! DONE
      CALL close_mesh(resformat,nres,ierr)
      CALL check_call(ierr,'GRETEL:CLOSEMESH:RES')
      CALL close_bnd(geoformat,ngeo,ierr)
      CALL check_call(ierr,'GRETEL:CLOSE_BND:GEO')
      CALL close_mesh(geoformat,ngeo,ierr)
      CALL check_call(ierr,'GRETEL:CLOSE_MESH:GEO')

      DEALLOCATE(textelu)
      WRITE(lu,*) 'END OF PROGRAM, ',ntimestep_res,' DATASETS FOUND'

      END SUBROUTINE gretel_autop