set_ps.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name: CHECK_DATA_09                                          C
!  Purpose: Check point source specifications.                         C
!                                                                      C
!  Author: J. Musser                                  Date: 10-JUN-13  C
!                                                                      C
!  Literature/Document References:                                     C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE set_ps

      use compar
      use exit, only: mfix_exit
      use functions
      use geometry
      use mpi_utility
      use param
      use param1, only: zero, small_number, undefined
      use physprop
      use ps
      use run

      implicit none

      INTEGER :: IJK, I, J, K, M, NN

      INTEGER PSV

      CHARACTER(LEN=64) :: eMsg

      INTEGER :: PS_SIZE

      DOUBLE PRECISION, allocatable :: lData_dp(:)
      DOUBLE PRECISION, allocatable :: gData_dp(:)

      LOGICAL, parameter :: dbg_PS = .false.

      if(.NOT.point_source) return

! DETERMINE WHICH BOUNDARY CONDITION INDICES HAVE VALUES
      l50: do psv = 1, dimension_ps

         IF(.NOT.ps_defined(psv)) cycle l50


! Calculate the velocity magnitude and normalize the axial components.
         CALL calc_ps_vel_mag(ps_vel_mag_g(psv), ps_u_g(psv),          &
            ps_v_g(psv), ps_w_g(psv))

         CALL calc_ps_cpxmflow(ps_cpxmflow_g(psv), ps_massflow_g(psv), &
            ps_t_g(psv), ps_x_g(psv,:), 0, c_pg0, dim_n_g, mw_g)

         do m=1, mmax
            CALL calc_ps_vel_mag(ps_vel_mag_s(psv,m), ps_u_s(psv,m),   &
               ps_v_s(psv,m), ps_w_s(psv,m))

            CALL calc_ps_cpxmflow(ps_cpxmflow_s(psv,m),                &
               ps_massflow_s(psv,m), ps_t_s(psv,m), ps_x_s(psv,m,:), m,&
               c_ps0(m), dim_n_s, mw_s(m,:))
         enddo


! Calculate the number of cells comprising the point source. This
! information is used to allocate some temp arrays.
!---------------------------------------------------------------------//
         ps_size = (ps_i_e(psv) - ps_i_w(psv) + 1) * &
                   (ps_j_n(psv) - ps_j_s(psv) + 1)
         if(do_k) ps_size = ps_size * (ps_k_t(psv) - ps_k_b(psv) + 1)

         if(ps_size < 1) then
             emsg = ''; write(emsg,"('Invalid PS size: ', I4)")ps_size
             goto 500
         endif


         allocate(ldata_dp( 0:numpes-1)); ldata_dp = zero
         allocate(gdata_dp( 0:numpes-1)); gdata_dp = zero


! Calculate the volume of the PointSource cells
!---------------------------------------------------------------------//
! Initialize the loop counter
         ps_volume(psv) = zero

         do k = ps_k_b(psv), ps_k_t(psv)
         do j = ps_j_s(psv), ps_j_n(psv)
         do i = ps_i_w(psv), ps_i_e(psv)

            if(is_on_mype_owns(i, j, k)) then
            IF (dead_cell_at(i,j,k)) cycle  ! skip dead cells
               ijk = funijk(i,j,k)

               if(fluid_at(ijk)) &
                  ldata_dp(mype) = ldata_dp(mype) + vol(ijk)
            endif

         enddo
         enddo
         enddo

! Each process (in DMP) only knows about the volume of the point source
! it sees. Invoke send/recv so that all process can calculate the total
! volume.
         CALL global_all_sum(ldata_dp, gdata_dp)

         ps_volume(psv) = sum(gdata_dp)
         if(ps_volume(psv) == zero) then
            emsg = 'No PS_VOLUME == ZERO'
            CALL debug_ps(psv, ps_size)
            goto 501
         endif

         if(allocated(ldata_dp)) deallocate(ldata_dp)
         if(allocated(gdata_dp)) deallocate(gdata_dp)


         IF(dbg_ps) CALL debug_ps(psv, ps_size)

      enddo l50

      return

  500 continue
      if(mype == pe_io) &
         write(*,"('PointSource setup Error: ',A)") trim(emsg)

      call mfix_exit(mype)

  501 continue
      write(*,"('PointSource setup Error: ',I3,2x,A)") mype, trim(emsg)

      call mfix_exit(mype)

      RETURN

      contains

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name: CALC_PS_CpxMFLOW                                       C
!  Purpose: Calculate velocity magnitude and normalize U/V/W.          C
!                                                                      C
!  Author: J. Musser                                  Date: 10-JUN-13  C
!                                                                      C
!  Literature/Document References:                                     C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE calc_ps_cpxmflow(CpxMFLOW, PS_MFLOW, PS_T, PS_X, lM, &
         cp0, ldim_n, lmw)

      use constant, only: gas_const_cal
      use read_thermochemical, only: calc_cpor

      DOUBLE PRECISION, intent(out) :: CpxMFLOW

      INTEGER, intent(in) :: lM
      INTEGER, intent(in) :: lDIM_N

      DOUBLE PRECISION, intent(in) :: Cp0
      DOUBLE PRECISION, intent(in) :: PS_MFLOW
      DOUBLE PRECISION, intent(in) :: PS_T          ! Temperature
      DOUBLE PRECISION, intent(in) :: PS_X(ldim_n)  ! Species Mass Frac

      DOUBLE PRECISION, intent(in) :: lMW(ldim_n)

      INTEGER :: IER

! If there is no mass flow for this phase, then there is no need to
! calculate a CPxMFLUX. Set it to zero and return.
      if(.NOT.energy_eq .OR. ps_mflow == zero) then
         cpxmflow = zero
         return
      endif

! Calculate the average specific heat.
      if(cp0 == undefined) then
         IF(.NOT.database_read) call read_database0()
         cpxmflow = zero
         do nn = 1, nmax(lm)
            cpxmflow = cpxmflow + ps_x(nn) * (gas_const_cal / lmw(nn)) * &
              calc_cpor(ps_t, lm, nn)
         enddo
      else
         cpxmflow = cp0
      endif

! Unit conversion (if needed)
      if(units == 'SI') cpxmflow = 4.183925d3*cpxmflow
      cpxmflow = cpxmflow*ps_mflow

      END SUBROUTINE calc_ps_cpxmflow

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name: CALC_PS_VEL_MAG                                        C
!  Purpose: Calculate velocity magnitude and normalize U/V/W.          C
!                                                                      C
!  Author: J. Musser                                  Date: 10-JUN-13  C
!                                                                      C
!  Literature/Document References:                                     C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE calc_ps_vel_mag(VEL_MAG, lU, lV, lW)

      DOUBLE PRECISION, intent(inout) :: VEL_MAG
      DOUBLE PRECISION, intent(inout) :: lU, lV, lW

! Normalize velocities:
      vel_mag = lu**2 + lv**2
      if(do_k) vel_mag = vel_mag + lw**2

      vel_mag = sqrt(vel_mag)

      if(vel_mag > small_number) then
         lu = lu/vel_mag
         lv = lv/vel_mag
         lw = lw/vel_mag
      else
         vel_mag = zero
         lu = zero
         lv = zero
         lw = zero
      endif


      END SUBROUTINE calc_ps_vel_mag


      END SUBROUTINE set_ps



!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name: DEBUG_PS                                               C
!                                                                      C
!  Purpose: Write out some information about that point source that    C
!  may be useful in debugging.                                         C
!                                                                      C
!  Author: J. Musser                                  Date: 24-JUN-13  C
!                                                                      C
!  Literature/Document References:                                     C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE debug_ps(lPSV, lPS_SIZE)

      use bc
      use compar
      use constant
      use cutcell
      use fldvar
      use geometry
      use ic
      use indices
      use mflux
      use parallel
      use param
      use param1
      use physprop
      use run
      use sendrecv
      use toleranc
      use usr
      use rxns
      use ps
      use mpi_utility
      use functions

      IMPLICIT NONE

! Index of PS source to debug.
      INTEGER, intent(in) :: lPSV
! Number of cells comprising the point source.
      INTEGER, intent(in) :: lPS_SIZE

      INTEGER :: IJK, I, J, K, M, NN

      INTEGER :: lc1

      INTEGER, allocatable :: lFlags_i(:,:)
      INTEGER, allocatable :: gFlags_i(:,:)

      if(mype == pe_io) then
         write(*,"(3/,3x,'Debug Point Source Index: ',I3)") lpsv
         write(*,"(/3x,'Size: ',I4)") lps_size
      endif

      allocate(lflags_i(lps_size,1:2) );   lflags_i = 0
      allocate(gflags_i(lps_size,1:2) );   gflags_i = 0

      lc1 = 0

      do k = ps_k_b(lpsv), ps_k_t(lpsv)
      do j = ps_j_s(lpsv), ps_j_n(lpsv)
      do i = ps_i_w(lpsv), ps_i_e(lpsv)

         lc1 = lc1 + 1
         if(is_on_mype_owns(i, j, k)) then
         IF (dead_cell_at(i,j,k)) cycle  ! skip dead cells
            ijk = funijk(i,j,k)
            if(fluid_at(ijk)) then
               lflags_i(lc1,1) = mype
               lflags_i(lc1,2) = flag(ijk)
            endif
         endif
      enddo
      enddo
      enddo

! Collect flag information on root.
      CALL global_sum(lflags_i, gflags_i)

! Write some information to the screen.
      if(mype == pe_io) then
         write(*,"(/5x,'Location:')")
         write(*,"( 5x,'X:',2(2x,g12.5),' :: ',2(2x,I4))")&
            ps_x_w(lpsv), ps_x_e(lpsv), ps_i_w(lpsv), ps_i_e(lpsv)
         write(*,"( 5x,'Y:',2(2x,g12.5),' :: ',2(2x,I4))")&
            ps_y_s(lpsv), ps_y_n(lpsv), ps_j_s(lpsv), ps_j_n(lpsv)
         if(do_k)write(*,"( 5x,'Z:',2(2x,g12.5),' :: ',2(2x,I4))")&
            ps_z_b(lpsv), ps_z_t(lpsv), ps_k_b(lpsv), ps_k_t(lpsv)

         write(*,"(/5x,'Volume: ',g12.5)") ps_volume(lpsv)


         if(ps_massflow_g(lpsv) > small_number) then
            write(*,"(//5x,'Point Source Gas Phase:')")
            write(*,"(7x,'Mass Flow Rate: ',g12.5)")ps_massflow_g(lpsv)
            write(*,"(7x,'Velocity Magnitude: ',g12.5)") ps_vel_mag_g(lpsv)
            write(*,"(7x,'Normal:')")
            write(*,"(9x,'x-Axis: ',g12.5)")ps_u_g(lpsv)
            write(*,"(9x,'y-Axis: ',g12.5)")ps_v_g(lpsv)
            if(do_k) write(*,"(9x,'z-Axis: ',g12.5)")ps_w_g(lpsv)
            if(energy_eq) &
               write(*,"(7x,'Temperature: ',g12.5)")ps_t_g(lpsv)
            if(species_eq(0)) then
               write(*,"(7x,'Species Composition:')")
               do nn=1,nmax(0)
                  write(*,"(9x,A,': ',g12.5)") &
                     trim(species_alias_g(nn)), ps_x_g(lpsv,nn)
               enddo
            endif
         else
            write(*,"(//5x,'No gas phase point source.')")
         endif


         do m=1,mmax
            if(ps_massflow_s(lpsv,m) > small_number) then
               write(*,"(//5x,'Point Source Solids Phase ',I1,':')")m
               write(*,"(7x,'Mass Flow Rate: ',g12.5)")ps_massflow_s(lpsv,m)
               write(*,"(7x,'Velocity Magnitude: ',g12.5)") ps_vel_mag_s(lpsv,m)
               write(*,"(7x,'Normal:')")
               write(*,"(9x,'x-Axis: ',g12.5)")ps_u_s(lpsv,m)
               write(*,"(9x,'y-Axis: ',g12.5)")ps_v_s(lpsv,m)
               if(do_k) write(*,"(9x,'z-Axis: ',g12.5)")ps_w_s(lpsv,m)
               if(energy_eq) &
                  write(*,"(7x,'Temperature: ',g12.5)")ps_t_s(lpsv,m)
               if(species_eq(m)) then
                  write(*,"(7x,'Species Composition:')")
                  do nn=1,nmax(m)
                     write(*,"(9x,A,': ',g12.5)") &
                        trim(species_alias_s(m,nn)), ps_x_s(lpsv,m,nn)
                  enddo
               endif
            else
               write(*,"(//5x,'No solids phase ',I1,' point source.')") m
            endif
         enddo

         write(*,"(//5x,'Point Source Cells:')")
         write(*,"(9x,'IJK',3(6x,A1),3x,'OWNS',3x,'FLAG')") 'I','J','K'

         lc1 = 0
         do k = ps_k_b(lpsv), ps_k_t(lpsv)
         do j = ps_j_s(lpsv), ps_j_n(lpsv)
         do i = ps_i_w(lpsv), ps_i_e(lpsv)
            IF (dead_cell_at(i,j,k)) cycle  ! skip dead cells
            lc1 = lc1 + 1
            write(*,"(4x,I8,5(3x,I4))") ijk, i, j, k,  gflags_i(lc1,:)
         enddo
         enddo
         enddo

      endif

      if(allocated(lflags_i)) deallocate(lflags_i)
      if(allocated(gflags_i)) deallocate(gflags_i)

      END SUBROUTINE debug_ps