write_des_data.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: WRITE_DES_DATA                                          !
!  Purpose: Writing DES output in Paraview format                      !
!                                                                      !
!                                                                      !
!  Author: Jay Boyalakuntla                           Date: 26-Jul-06  !
!  Reviewer: Sreekanth Pannala                        Date: 31-Oct-06  !
!                                                                      !
!  Reviewer: J. Musser                                Date: 20-Apr-10  !
!  Comments: Split original subroutine into one for ParaView *.vtp     !
!  files, and a second for TECPLOT files *.dat.                        !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!

      SUBROUTINE write_des_data

!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE param
      USE param1
      USE parallel
      USE fldvar
      USE discretelement
      USE run
      USE geometry
      USE physprop
      USE sendrecv
      USE des_bc

      use error_manager

      IMPLICIT NONE
!-----------------------------------------------
! Local Variables
!-----------------------------------------------

!-----------------------------------------------
! Functions
!-----------------------------------------------

!-----------------------------------------------

      IF (trim(des_output_type) .EQ. 'TECPLOT') THEN
         CALL write_des_tecplot
      ELSE
         CALL write_des_vtp
      ENDIF

! Invoke at own risk

! Granular temperature subroutine should be called/calculated when
! writing DES data

      IF (des_calc_bedheight) THEN
         CALL calc_des_bedheight
         CALL write_des_bedheight
      ENDIF

      IF (.false.) CALL des_granular_temperature
      IF (.false.) CALL write_des_theta

      RETURN
      END SUBROUTINE write_des_data
!-----------------------------------------------

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: WRITE_DES_VTP                                           !
!  Purpose: Writing DES output in Paraview format                      !
!                                                                      !
!                                                                      !
!  Reviewer: Rahul Garg                               Dare: 01-Aug-07  !
!  Comments: Added one more output file containing averaged bed height !
!                                                                      !
!  Revision : For parallel runs added cumulative and parallel IO       !
!  Author   : Pradeep G.                                               !
!                                                                      !
!  NOTE: If the system starts with zero particles, ParaView may have   !
!  trouble showing the results. To view the results in the current     !
!  version of ParaView, Version 3.6.1:                                 !
!    i - load the *.vtp files                                          !
!   ii - filter with glyph (Filters > Common > Glyph)                  !
!        a - change glyph to sphere                                    !
!        b - change scale mode to scalar                               !
!        c - check the "Edit" Set Scale Factor Box                     !
!        d - change the value to 1.0                                   !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!

      SUBROUTINE write_des_vtp

      use des_rxns, only: des_x_s
      use des_thermo, only: des_t_s
      use discretelement, only: des_pos_new, des_vel_new, des_usr_var
      use discretelement, only: des_radius
      use discretelement, only: des_usr_var, des_usr_var_size
      use discretelement, only: s_time
      use discretelement, only: use_cohesion, postcohesive
      use error_manager
      use mfix_pic, only: des_stat_wt, mppic
      use param
      use param1, only: zero
      use run, only: any_species_eq
      use run, only: energy_eq
      use vtp

      IMPLICIT NONE

      DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: TMP_PAR
      CHARACTER(len=10) :: lNoP
      CHARACTER(len=24) :: sTIMEc
      INTEGER :: NN

      stimec=''; WRITE(stimec,"(ES24.16)") s_time

! This routine opens the VTP file and calculates send/recv information.
! It returns back the number of points as a string.
      CALL vtp_open_file(lnop)

! Standard VTP header information:
!----------------------------------------------------------------------/
      CALL vtp_write_element('<?xml version="1.0"?>')
      CALL vtp_write_element('<!-- Time ='//stimec//'s -->')
      CALL vtp_write_element('<VTKFile type="PolyData" version="0.1" &
         &byte_order="LittleEndian">')
      CALL vtp_write_element('<PolyData>')

      CALL vtp_write_element('<Piece NumberOfPoints="'//lnop//'" &
         &NumberOfVerts="0" NumberOfLines="0" NumberOfStrips="0" &
         &NumberOfPolys="0">')

! Points are the particle identified by position:
!----------------------------------------------------------------------/
      CALL vtp_write_element('<Points>')
      CALL vtp_write_data('Position', des_pos_new)
      CALL vtp_write_element('</Points>')

! PointData are individual particle properties:
!----------------------------------------------------------------------/
      ALLOCATE(tmp_par(SIZE(des_radius)))

      CALL vtp_write_element('<PointData Scalars="Diameter" &
         &Vectors="Velocity">')

      CALL get_diameter(tmp_par)
      CALL vtp_write_data('Diameter', tmp_par)

      CALL vtp_write_data('Velocity', des_vel_new)

      CALL vtp_write_data('ID', iglobal_id)

      DO nn=1, des_usr_var_size
         CALL vtp_write_data(trim(ivar('USR_VAR',nn)),des_usr_var(nn,:))
      ENDDO

      IF(particle_orientation) &
         CALL vtp_write_data('Orientation', orientation)

      IF(mppic) THEN
         tmp_par = 1.0d0 - epg_p
         CALL vtp_write_data('EPs', tmp_par)
      ENDIF

      IF(energy_eq) &
         CALL vtp_write_data('Temperature', des_t_s)

      IF(any_species_eq) THEN
         DO nn=1, dimension_n_s
            CALL vtp_write_data(trim(ivar('X_s',nn)), des_x_s(:,nn))
         ENDDO
      ENDIF

      IF(use_cohesion) &
         CALL vtp_write_data('CohesiveForce', postcohesive)

      CALL vtp_write_element('</PointData>')
      DEALLOCATE(tmp_par)

! Open/Close the unused VTP tags.
!----------------------------------------------------------------------/
      CALL vtp_write_element('<CellData></CellData>')
      CALL vtp_write_element('<Verts></Verts>')
      CALL vtp_write_element('<Lines></Lines>')
      CALL vtp_write_element('<Strips></Strips>')
      CALL vtp_write_element('<Polys></Polys>')

! Close all the opened tags:
!----------------------------------------------------------------------/
      CALL vtp_write_element('</Piece>')
      CALL vtp_write_element('</PolyData>')
      CALL vtp_write_element('</VTKFile>')

      CALL vtp_close_file

! Add the new VTP file to the PVD file for time association.
      CALL add_vtp_to_pvd

      RETURN

      contains

!......................................................................!
! Purpose: Return the diamerter of DES particles.                      !
!......................................................................!
      SUBROUTINE get_diameter(OUT_VAR)

      IMPLICIT NONE

      DOUBLE PRECISION, INTENT(OUT) :: OUT_VAR(:)
      INTEGER :: LC

! Return the effective diamter of the parcle.
      IF(mppic) THEN
         DO lc=1,size(out_var)
            IF(des_stat_wt(lc) > zero) THEN
               out_var(lc) = 2.0d0*des_radius(lc)*&
                  (des_stat_wt(lc)**(1./3.))
            ELSE
               out_var(lc) = 2.0d0*des_radius(lc)
            ENDIF
         ENDDO
      ELSE
         out_var = 2.0d0*des_radius
      ENDIF

      RETURN
      END SUBROUTINE get_diameter

      END SUBROUTINE write_des_vtp

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
!
!  Module name: WRITE_DES_TECPLOT
!  Purpose: Writing DES output in TECPLOT format
!
!  Revision: For parallel runs added distributed and single IO
!  Comment: In earlier version the time instances are keep appended to
!           the tecplot file. This will make tecplot file so large for
!           large simulations. Hence seperate files are written for
!           each instances
!  Author : Pradeep G.
!
!
!
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      SUBROUTINE write_des_tecplot

      USE param
      USE param1
      USE parallel
      USE fldvar
      USE discretelement
      USE run
      USE geometry
      USE physprop
      USE sendrecv
      USE des_bc
      USE compar
      USE cdist
      USE desmpi
      USE mpi_comm_des
      USE mpi_utility
      USE functions
      IMPLICIT NONE

!-----------------------------------------------
! Local Variables
!-----------------------------------------------
! file units for
! Pradeep remove parameter from following variables
      INTEGER:: DES_DATA,DES_EX,DES_EPS

! output file for basic DES variables including: position, velocity,
! radius, density, mark (flag)
      CHARACTER(LEN=50)     :: FNAME_DATA

! output file for extra DES variables including:
! solids time (S_TIME), maximum neighbor count, maximum overlap
! granular energy and granular temperature
      CHARACTER(LEN=50)     :: FNAME_EXTRA

! output file for axial solids volume fraction and granular temp
      CHARACTER(LEN=50)     :: FNAME_EPS

! dummy indices
      INTEGER L, K

! index to track accounted for particles
      INTEGER PC

! Variables related to gathering info at PE_IO
      integer llocalcnt,lglocnt,lgathercnts(0:numpes-1),lproc,ltotvar,lcount
      double precision, dimension(:,:), allocatable :: ltemp_array

      INTEGER :: wDIMN

! Set output dimnensions
      wdimn = merge(2,3,no_k)
! set the total variable based on dimension
      ltotvar = merge(8,9,no_k)

! set the file name and unit number and open file
      des_data = 2000
      des_ex = 2100
      des_eps = 2200
      if (bdist_io) then
         write(fname_data,'(A,"_DES_DATA",I4.4,"_",I4.4,".dat")') trim(run_name),tecplot_findex,mype
         write(fname_extra,'(A,"_DES_EXTRA",I4.4,"_",I4.4,".dat")') trim(run_name),tecplot_findex,mype
         write(fname_eps,'(A,"_DES_EPS",I4.4,"_",I4.4,".dat")') trim(run_name),tecplot_findex,mype
         open(convert='big_endian',unit=des_data,file=fname_data,status='new',err=999)
         open(convert='big_endian',unit=des_ex,file=fname_extra,status='new',err=999)
         open(convert='big_endian',unit=des_eps,file=fname_eps,status='new',err=999)
      else
         if(mype.eq.pe_io) then
            write(fname_data,'(A,"_DES_DATA_",I4.4,".dat")') trim(run_name),tecplot_findex
            write(fname_extra,'(A,"_DES_EXTRA_",I4.4,".dat")') trim(run_name),tecplot_findex
            write(fname_eps,'(A,"_DES_EPS_",I4.4,".dat")') trim(run_name),tecplot_findex
            open(convert='big_endian',unit=des_data,file=fname_data,status='new',err=999)
            open(convert='big_endian',unit=des_ex,file=fname_extra,status='new',err=999)
            open(convert='big_endian',unit=des_eps,file=fname_eps,status='new',err=999)
         end if
      end if
      tecplot_findex = tecplot_findex + 1

! write header
      if (bdist_io .or. mype .eq. pe_io) then
         if(do_k) then
            write (des_data,'(A)') &
            'variables = "x" "y" "z" "vx" "vy" "vz" "rad" "den" "mark"'
         else
            write (des_data, '(A)') &
            'variables = "x" "y" "vx" "vy" "omega" "rad" "den" "mark"'
         endif
         write (des_data, "(A,F15.7,A)")'zone T ="', s_time, '"'
      end if

! write data in ditributed mode or as a single file
      if(bdist_io) then
         pc = 1
         do l = 1,max_pip
            if(pc.gt.pip) exit
            if(is_nonexistent(l)) cycle
            pc = pc+1
            if(is_ghost(l) .or. is_entering_ghost(l) .or. is_exiting_ghost(l)) cycle
            if(do_k) then
               write (des_data, '(8(2x,es12.5))')&
                  (des_pos_new(l,k),k=1,wdimn),(des_vel_new(l,k),k=1,wdimn), &
                   des_radius(l),ro_sol(l)
            else
               write (des_data, '(7(2x,es12.5))')&
                  (des_pos_new(l,k),k=1,wdimn), (des_vel_new(l,k),k=1,wdimn), &
                  omega_new(l,1), des_radius(l), ro_sol(l)
            endif
        end do

      else ! if bdist_io
! set parameters required for gathering info at PEIO and write as single file
         lglocnt = 10
         llocalcnt = pip - ighost_cnt
         call global_sum(llocalcnt,lglocnt)
         allocate (dprocbuf(llocalcnt),drootbuf(lglocnt),iprocbuf(llocalcnt),irootbuf(lglocnt))
         allocate (ltemp_array(lglocnt,ltotvar))
         igath_sendcnt = llocalcnt
         lgathercnts = 0
         lgathercnts(mype) = llocalcnt
         call global_sum(lgathercnts,igathercnts)
         idispls(0) = 0
         do lproc = 1,numpes-1
            idispls(lproc) = idispls(lproc-1) + igathercnts(lproc-1)
         end do

! gather information from all processor
         lcount = 1
         do k = 1,wdimn
            call des_gather(des_pos_new(:,k))
            ltemp_array(:,lcount) = drootbuf(:); lcount=lcount+1
         end do
         do k = 1,wdimn
            call des_gather(des_vel_new(:,k))
            ltemp_array(:,lcount) = drootbuf(:); lcount=lcount+1
         end do
         if(no_k) then
            call des_gather(omega_new(:,1))
            ltemp_array(:,lcount) = drootbuf(:); lcount=lcount+1
         end if
         call des_gather(des_radius)
         ltemp_array(:,lcount) = drootbuf(:); lcount=lcount+1
         call des_gather(ro_sol)
         ltemp_array(:,lcount) = drootbuf(:); lcount=lcount+1

! write the data into file
         if (mype.eq.pe_io) then
            if(do_k) then
               do l =1,lglocnt
                  write (des_data,'(8(2x,es12.5),I5)') (ltemp_array(l,k),k=1,8),int(ltemp_array(l,9))
               end do
            else
               do l =1,lglocnt
                  write (des_data,'(7(2x,es12.5),I5)') (ltemp_array(l,k),k=1,7),int(ltemp_array(l,8))
               end do
            end if
         end if
         deallocate (dprocbuf,drootbuf,iprocbuf,irootbuf,ltemp_array)
      end if
! close the files
      if (mype.eq.pe_io .or. bdist_io) then
         close(des_data)
         close(des_ex)
         close(des_eps)
      end if
      return

  999 write(*,"(/1x,70('*'),//,a,/,a,/1x,70('*'))")&
         ' From: write_des_tecplot ',&
         ' message: error opening des tecplot file. terminating run.'


      END SUBROUTINE write_des_tecplot
!-----------------------------------------------

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
!
!  Module name: WRITE_DES_BEDHEIGHT
!  Purpose: Writing DES output on bed height.

!  WARNING: This code is out-of-date and should be modified for consistency
!  with current DEM version.  Also this routine will be fairly specific
!  to a user needs and should probably be tailored as such

!
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      SUBROUTINE write_des_bedheight

      USE param
      USE param1
      USE parallel
      USE fldvar
      USE discretelement
      USE run
      USE geometry
      USE physprop
      USE sendrecv
      USE des_bc
      IMPLICIT NONE

!-----------------------------------------------
! Local Variables
!-----------------------------------------------
! logical that identifies that the data file has been created
! and is already opened (initial checks/writes)
      LOGICAL, SAVE :: FIRST_PASS = .true.
! logical used for testing is the data files already exists
      LOGICAL :: F_EXISTS
! output file for the bed height data
      CHARACTER(LEN=50)     :: FNAME_BH
! file unit for the bed height data
      INTEGER, PARAMETER :: BH_UNIT = 2010
! dummy index values
      INTEGER I, M
! variables for bed height calculation
      INTEGER, SAVE :: tcount = 1
      DOUBLE PRECISION :: height_avg, height_rms
      DOUBLE PRECISION, PARAMETER :: tmin = 5.d0
      DOUBLE PRECISION, DIMENSION(5000), SAVE :: bed_height_time

!-----------------------------------------------
! Functions
!-----------------------------------------------
!-----------------------------------------------

! after tmin start storing bed height. after enough measurements
! have been taken (i.e. tcount > 20) start to calculate a running
! average bed height and running rms bed height for solids phase 1 only
      height_avg = zero
      height_rms = zero

      if(time.gt.tmin) then
         if(tcount.le.5000)  then
            bed_height_time(tcount) = bed_height(1)
            !dt_time(tcount) = DT
            tcount = tcount + 1

            if(tcount.gt.20)  then
               do i = 1, tcount-1,1
                  height_avg = height_avg + bed_height_time(i)!*dt_time(i)
               enddo
               height_avg = height_avg/(tcount-1)
               do i = 1, tcount-1,1
                  height_rms = height_rms + ((bed_height_time(i)&
                       &-height_avg)**2)!*dt_time(i)
               enddo

               height_rms = sqrt(height_rms/(tcount-1))
            endif
         endif
      endif

      fname_bh = trim(run_name)//'_DES_BEDHEIGHT.dat'
      IF(first_pass) THEN
         f_exists = .false.
         INQUIRE(file=fname_bh,exist=f_exists)
! If the file does not exist, then create it with the necessary
! header information.
         IF (.NOT.f_exists) THEN
            OPEN(convert='BIG_ENDIAN',unit=bh_unit,file=fname_bh,&
               form="formatted",status="new")
         ELSE
! To prevent overwriting existing files accidently, exit if the file
! exists and this is a NEW run.
            IF(run_type .EQ. 'NEW') THEN
               WRITE(*,1000)
               WRITE(unit_log, 1000)
               CALL mfix_exit(mype)
            ELSE
! Open the file for appending of new data (RESTART_1 Case)
               OPEN(convert='BIG_ENDIAN',unit=bh_unit,file=fname_bh,&
                    position="append")
            ENDIF
         ENDIF
         first_pass = .false.
      ELSE
! Open the file and mark for appending
         OPEN(convert='BIG_ENDIAN',unit=bh_unit,file=fname_bh,&
              position="append")
      ENDIF

      WRITE(bh_unit, '(10(2X,E20.12))') s_time, &
         (bed_height(m), m=mmax+1,mmax+des_mmax), height_avg, height_rms
! Close the file and keep
      CLOSE(bh_unit, status="KEEP")

      RETURN

 1000 FORMAT(/1x,70('*')//, ' From: WRITE_DES_BEDHEIGHT',/,&
         ' Message: bed_height.dat already exists in the run',&
         ' directory.',/10x, 'Run terminated to prevent',&
         ' accidental overwriting of files.',/1x,70('*')/)

      END SUBROUTINE write_des_bedheight


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
!
!  Module name: WRITE_DES_THETA
!  Purpose: The following code writes out theta_m for discrete
!  particles to a file for each ijk cell in the system each time
!  des_granular_temperature is called.
!
!
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      SUBROUTINE write_des_theta

      USE param
      USE param1
      USE parallel
      USE fldvar
      USE discretelement
      USE run
      USE geometry
      USE physprop
      USE sendrecv
      USE des_bc
      USE functions
      IMPLICIT NONE

!-----------------------------------------------
! Local Variables
!-----------------------------------------------
! indices
      INTEGER I, J, K, IJK
!
      INTEGER M, NP
! logical that identifies that the data file has been created
! and is already opened (initial checks/writes)
      LOGICAL, SAVE :: FIRST_PASS = .true.
! logical used for testing is the data file already exists
      LOGICAL :: F_EXISTS
! file unit for the granular temperature data
      INTEGER, PARAMETER :: GT_UNIT = 2020
! output file for the granular temperature data
      CHARACTER(LEN=50)  :: FNAME_GT
!-----------------------------------------------

      fname_gt = trim(run_name)//'_DES_THETA.dat'
      IF (first_pass) THEN
         f_exists = .false.
         INQUIRE(file=fname_gt,exist=f_exists)

         IF (.NOT.f_exists) THEN
! If the file does not exist, then create it with the necessary
! header information.
            OPEN(convert='BIG_ENDIAN',unit=gt_unit,file=fname_gt,&
                 status='NEW')
         ELSE
            IF(run_type .EQ. 'NEW') THEN
! If the run is new and the GT file already exists replace it with a
! new file.
!               OPEN(CONVERT='BIG_ENDIAN',UNIT=GT_UNIT,FILE=FNAME_GT,STATUS='REPLACE')
! Prevent overwriting an existing file by exiting if the file exists
! and this is a NEW run.
               WRITE(*,1001) fname_gt
               WRITE(unit_log,1001) fname_gt
               CALL mfix_exit(mype)
            ELSE
! Open the file for appending of new data (RESTART_1 Case)
               OPEN(convert='BIG_ENDIAN',unit=gt_unit, file=fname_gt,&
                    position='APPEND')
            ENDIF
         ENDIF
         first_pass =  .false.
      ELSE
! Open file and mark for appending
         OPEN(convert='BIG_ENDIAN',unit=gt_unit,file=fname_gt,&
              position='APPEND')
      ENDIF   ! endif (first_pass)

      WRITE(gt_unit,*) ''
      WRITE(gt_unit,'(A6,ES24.16)') 'Time=', s_time
      WRITE(gt_unit,'(A6,2X,3(A6,2X),A8)',advance="NO") 'IJK', &
         'I', 'J', 'K', 'NP'
      DO m = mmax+1,des_mmax+mmax
         WRITE(gt_unit,'(7X,A6,I1)',advance="NO") 'THETA_',m
      ENDDO
      WRITE(gt_unit,*) ''
      DO ijk = ijkstart3, ijkend3
         IF(fluid_at(ijk)) THEN
            i = i_of(ijk)
            j = j_of(ijk)
            k = k_of(ijk)
            np = pinc(ijk)
            WRITE(gt_unit,'(I6,2X,3(I6,2X),I8,(2X,ES15.5))') &
               ijk, i, j, k, np, (theta_m(ijk,m), m = mmax+1,des_mmax+mmax)
         ENDIF
      ENDDO

! Close the file and keep
      CLOSE(gt_unit, status='KEEP')

      RETURN

 1001 FORMAT(/1x,70('*')//, ' From: WRITE_DES_THETA',/,&
         ' Message: ', a, ' already exists in the run',/10x,&
         'directory. Run terminated to prevent accidental overwriting',&
         /10x,'of files.',/1x,70('*')/)

      END SUBROUTINE write_des_theta
!-----------------------------------------------