solve_granular_energy.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: SOLVE_GRANULAR_ENERGY                                   C
!  Purpose: Solve granular energy equations in matrix equation         C
!     form Ax=b. The center coefficient (ap) and source vector (b)     C
!     are negative.  The off-diagonal coefficients are positive.       C
!                                                                      C
!  Author: K. Agrawal                                 Date:            C
!                                                                      C
!  Literature/Document References:                                     C
!                                                                      C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE solve_granular_energy(IER)

! Modules
!---------------------------------------------------------------------//
      use ambm, only: a_m, b_m, lock_ambm, unlock_ambm
      use bc, only: bc_theta_m, bc_thetaw_m, bc_hw_theta_m, bc_c_theta_m
      use compar, only: ijkstart3, ijkend3, mype, numpes
      use constant, only: to_si, pi
      use exit, only: mfix_exit
      use fldvar, only: ep_g, theta_m, theta_mo
      use fldvar, only: ep_s, u_s, v_s, w_s, rop_so
      use fldvar, only: ro_s, d_p
      use functions, only: fluid_at, wall_at, zmax
      use geometry, only: ijkmax2, vol
      use leqsol, only: leq_it, leq_method, leq_sweep, leq_pc, leq_tol
      use mflux, only: flux_se, flux_sn, flux_st
      use mflux, only: flux_sse, flux_ssn, flux_sst
      use mflux, only: flux_ne, flux_nn, flux_nt
      use mms, only: use_mms, mms_theta_m_src
      use mpi_utility, only: global_all_sum
      use param, only: dimension_3, dimension_m
      use param1, only: zero, one, dimension_lm
      use physprop, only: kth_s
      use physprop, only: smax, mmax
      use residual, only: resid, max_resid, ijk_resid
      use residual, only: num_resid, den_resid
      use residual, only: resid_th
      use run, only: discretize, odt, added_mass, m_am, schaeffer
      use run, only: kt_type_enum, lun_1984, ahmadi_1995, simonin_1996
      use run, only: kt_type, gd_1999, gtsh_2012, ghd_2007, ia_2005
      use rxns, only: sum_r_s
      use toleranc, only: zero_ep_s 
      use ur_facs, only: ur_fac
      use usr_src, only: call_usr_source, calc_usr_source
      use usr_src, only: gran_energy
      use visc_s, only: ep_g_blend_start
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! Error index
      INTEGER :: IER

! Local variables
!---------------------------------------------------------------------//
! phase index
      INTEGER :: M, L
! Cp * Flux
      DOUBLE PRECISION CpxFlux_E(dimension_3), CpxFlux_N(dimension_3), CpxFlux_T(dimension_3)
! previous time step term
      DOUBLE PRECISION :: apo
! source terms which appear appear in the center coefficient (lhs) and
! the rhs vector
      DOUBLE PRECISION :: sourcelhs, sourcerhs
! Indices
      INTEGER :: IJK
! linear equation solver method and iterations
      INTEGER :: LEQM, LEQI
! particle mass and diameter
      DOUBLE PRECISION :: M_PM,D_PM
! total solids volume fraction
      DOUBLE PRECISION :: TOT_EPS(dimension_3)
! net production of all solids phase
      DOUBLE PRECISION :: TOT_SUM_RS(dimension_3)
! old value of total solids number density
      DOUBLE PRECISION :: TOT_NO(dimension_3)
! small value of theta_m, 1 cm2/s2 = 1e-4 m2/s2
      DOUBLE PRECISION :: smallTheta
! temporary variable for volume x interphase transfer
! coefficient between solids phases
      DOUBLE PRECISION :: VxTC_ss(dimension_3,dimension_lm)

! Arrays for storing errors:
! 140 - Linear Eq diverged
! 141 - Unphysical values
! 14x - Unclassified
      INTEGER :: Err_l(0:numpes-1)  ! local
      INTEGER :: Err_g(0:numpes-1)  ! global

! temporary use of global arrays:
! array1 (locally s_p)
! source vector: coefficient of dependent variable
! becomes part of a_m matrix; must be positive
      DOUBLE PRECISION :: S_P(dimension_3)
! array2 (locally s_c)
! source vector: constant part becomes part of b_m vector
      DOUBLE PRECISION :: S_C(dimension_3)
! array3 (locally eps)
      DOUBLE PRECISION :: eps(dimension_3)

! Septadiagonal matrix A_m, vector b_m
!      DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
!      DOUBLE PRECISION B_m(DIMENSION_3, 0:DIMENSION_M)
!---------------------------------------------------------------------//

      call lock_ambm       ! locks arrays a_m and b_m

! Initialize error flags.
      err_l = 0

      smalltheta = (to_si)**4 * zero_ep_s

      DO m = 0, mmax
         CALL init_ab_m (a_m, b_m, ijkmax2, m)
      ENDDO

      SELECT CASE (kt_type_enum)
        CASE(lun_1984, ahmadi_1995, simonin_1996, gd_1999, gtsh_2012)
! ---------------------------------------------------------------->>>
          DO m = 1, smax

            DO ijk = ijkstart3, ijkend3

               IF(.NOT.added_mass .OR. m /= m_am) THEN
                  cpxflux_e(ijk) = 1.5d0 * flux_se(ijk,m)
                  cpxflux_n(ijk) = 1.5d0 * flux_sn(ijk,m)
                  cpxflux_t(ijk) = 1.5d0 * flux_st(ijk,m)
               ELSE
                  cpxflux_e(ijk) = 1.5d0 * flux_sse(ijk)
                  cpxflux_n(ijk) = 1.5d0 * flux_ssn(ijk)
                  cpxflux_t(ijk) = 1.5d0 * flux_sst(ijk)
               ENDIF

               IF (fluid_at(ijk)) THEN
! calculate the source terms to be used in the a matrix and b vector
                  IF (kt_type_enum == gd_1999 .OR. &
                      kt_type_enum == gtsh_2012) THEN
                     CALL source_granular_energy_gd(sourcelhs, &
                        sourcerhs, ijk, m)
                  ELSE
                     CALL source_granular_energy(sourcelhs, &
                        sourcerhs, ijk, m)
                  ENDIF
                  apo = 1.5d0*rop_so(ijk,m)*vol(ijk)*odt
                  s_p(ijk) = apo + sourcelhs + 1.5d0* &
                     zmax(sum_r_s(ijk,m)) * vol(ijk)
                  s_c(ijk) = apo*theta_mo(ijk,m) + sourcerhs + 1.5d0* &
                      theta_m(ijk,m)*zmax((-sum_r_s(ijk,m))) * vol(ijk)
                  eps(ijk) = ep_s(ijk,m)
! MMS Source term.
                  IF(use_mms) s_c(ijk) = s_c(ijk) + &
                     mms_theta_m_src(ijk)*vol(ijk)
               ELSE
                  eps(ijk) = zero
                  s_p(ijk) = zero
                  s_c(ijk) = zero
                  IF(use_mms) eps(ijk) = ep_s(ijk,m)
               ENDIF
            ENDDO   ! end do loop (ijk=ijkstart3,ijkend3)

! calculate the convection-diffusion terms
            CALL conv_dif_phi (theta_m(1,m), kth_s(1,m), discretize(8), &
               u_s(1,m), v_s(1,m), w_s(1,m), &
               cpxflux_e, cpxflux_n, cpxflux_t, m, a_m, b_m)

! calculate standard bc
            CALL bc_phi (theta_m(1,m), bc_theta_m(1,m), bc_thetaw_m(1,m), &
               bc_hw_theta_m(1,m), bc_c_theta_m(1,m), m, a_m, b_m)

! override bc settings if Johnson-Jackson bcs are specified
            CALL bc_theta (m, a_m, b_m)

! set the source terms in a and b matrix form
            CALL source_phi (s_p, s_c, eps, theta_m(1,m), m, a_m, b_m)

! usr source
            IF (call_usr_source(8)) CALL calc_usr_source(gran_energy,&
                                  a_m, b_m, lm=m)

! Adjusting the values of theta_m to zero when Ep_g < EP_star
! (Shaeffer, 1987). This is done here instead of calc_mu_s to
! avoid convergence problems. (sof)
            IF (schaeffer) THEN
               DO ijk = ijkstart3, ijkend3
                  IF (fluid_at(ijk) .AND. &
                      ep_g(ijk) .LT. ep_g_blend_start(ijk)) THEN

                     a_m(ijk,1,m) = zero
                     a_m(ijk,-1,m) = zero
                     a_m(ijk,2,m) = zero
                     a_m(ijk,-2,m) = zero
                     a_m(ijk,3,m) = zero
                     a_m(ijk,-3,m) = zero
                     a_m(ijk,0,m) = -one
                     b_m(ijk,m) = -smalltheta
                  ENDIF
               ENDDO
            ENDIF

            CALL calc_resid_s (theta_m(1,m), a_m, b_m, m, &
               num_resid(resid_th,m), &
               den_resid(resid_th,m), resid(resid_th,m), &
               max_resid(resid_th,m), ijk_resid(resid_th,m), zero)

            CALL under_relax_s (theta_m(1,m), a_m, b_m, m, ur_fac(8))

!            call check_ab_m(a_m, b_m, m, .true., ier)
!            write(*,*) resid(resid_th, m), max_resid(resid_th, m),&
!               ijk_resid(resid_th, m)
!            call write_ab_m(a_m, b_m, ijkmax2, m, ier)

!            call test_lin_eq(ijkmax2, ijmax2, imax2, a_m(1, -3, M), &
!               1, DO_K, ier)

            CALL adjust_leq (resid(resid_th,m), leq_it(8), leq_method(8),&
               leqi, leqm)

            CALL solve_lin_eq ('Theta_m', 8, theta_m(1,m), a_m, b_m, m,&
               leqi, leqm, leq_sweep(8), leq_tol(8),  leq_pc(8), ier)
!            call out_array(Theta_m(1,m), 'Theta_m')

! Check for linear solver divergence.
            IF(ier == -2) err_l(mype) = 140

! Remove very small negative values of theta caused by leq solvers
            CALL adjust_theta (m, ier)
! large negative granular temp -> divergence
            IF (ier /= 0) err_l(mype) = 141

          ENDDO   ! end do loop (m=1,smax)
! end case(lun_1984, simonin_1996, ahmadi_1995, gd_1999, gtsh_2012)
! ----------------------------------------------------------------<<<

        CASE(ia_2005)
! ---------------------------------------------------------------->>>
          DO m = 1, smax
            DO ijk = ijkstart3, ijkend3

! Skip walls where some values are undefined.
               IF(wall_at(ijk)) cycle

               d_pm = d_p(ijk,m)
               m_pm = (pi/6.d0)*(d_pm**3)*ro_s(ijk,m)

! In Iddir & Arastoopour (2005) the granular temperature includes
! mass of the particle in the definition.
               IF(.NOT.added_mass .OR. m /= m_am) THEN
                  cpxflux_e(ijk) = (1.5d0/m_pm) * flux_se(ijk,m)
                  cpxflux_n(ijk) = (1.5d0/m_pm) * flux_sn(ijk,m)
                  cpxflux_t(ijk) = (1.5d0/m_pm) * flux_st(ijk,m)
               ELSE ! in case added mass is used.
                  cpxflux_e(ijk) = (1.5d0/m_pm) * flux_sse(ijk)
                  cpxflux_n(ijk) = (1.5d0/m_pm) * flux_ssn(ijk)
                  cpxflux_t(ijk) = (1.5d0/m_pm) * flux_sst(ijk)
               ENDIF

               IF (fluid_at(ijk)) THEN
! calculate the source terms to be used in the a matrix and b vector
                  CALL source_granular_energy_ia(sourcelhs, &
                     sourcerhs, ijk, m)
                  apo = (1.5d0/m_pm)*rop_so(ijk,m)*vol(ijk)*odt
                  s_p(ijk) = apo + sourcelhs + (1.5d0/m_pm)*&
                     zmax(sum_r_s(ijk,m)) * vol(ijk)
                  s_c(ijk) = apo*theta_mo(ijk,m) + sourcerhs + &
                     (1.50d0/m_pm)*theta_m(ijk,m)*&
                     zmax((-sum_r_s(ijk,m))) * vol(ijk)
                  eps(ijk) = ep_s(ijk,m)
               ELSE
                  eps(ijk) = zero
                  s_p(ijk) = zero
                  s_c(ijk) = zero
               ENDIF
            ENDDO    ! end do loop (ijk=ijkstart3,ijkend3)

! calculate the convection-diffusion terms
            CALL conv_dif_phi (theta_m(1,m), kth_s(1,m), discretize(8),&
               u_s(1,m), v_s(1,m), w_s(1,m), &
               cpxflux_e, cpxflux_n, cpxflux_t, m, a_m, b_m)

! calculate standard bc
            CALL bc_phi (theta_m(1,m), bc_theta_m(1,m), bc_thetaw_m(1,m),&
               bc_hw_theta_m(1,m),bc_c_theta_m(1,m), m, a_m, b_m)

! override bc settings if Johnson-Jackson bcs are specified
            CALL bc_theta (m, a_m, b_m)

! set the source terms in a and b matrix form
            CALL source_phi (s_p, s_c, eps, theta_m(1,m), m, a_m, b_m)

! usr source
            IF (call_usr_source(8)) CALL calc_usr_source(gran_energy,&
                                  a_m, b_m, lm=m)
          ENDDO   ! end do loop (m = 1, smax)

! use partial elimination on collisional dissipation term: SUM(Nip)
!     SUM( ED_s_ip* (Theta_p-Theta_i))
          IF (smax > 1) THEN
            CALL calc_vtc_ss (vxtc_ss)
            CALL partial_elim_ia (theta_m, vxtc_ss, a_m, b_m)
          ENDIF

! Adjusting the values of theta_m to zero when Ep_g < EP_star
! (Shaeffer, 1987). This is done here instead of calc_mu_s to
! avoid convergence problems. (sof)
          IF (schaeffer) THEN
            DO m = 1, smax
               DO ijk = ijkstart3, ijkend3
                  IF (fluid_at(ijk) .AND. &
                      ep_g(ijk) .LT. ep_g_blend_start(ijk)) THEN

                     d_pm = d_p(ijk,m)
                     m_pm = (pi/6.d0)*(d_pm**3)*ro_s(ijk,m)
                     a_m(ijk,1,m) = zero
                     a_m(ijk,-1,m) = zero
                     a_m(ijk,2,m) = zero
                     a_m(ijk,-2,m) = zero
                     a_m(ijk,3,m) = zero
                     a_m(ijk,-3,m) = zero
                     a_m(ijk,0,m) = -one
! In Iddir & Arastoopour (2005) the granular temperature includes
! mass of the particle in the definition. Systems with small
! particles can give rise to smaller temperatures than the standard
! zero_ep_s
                     b_m(ijk,m) = -smalltheta*m_pm
                  ENDIF
               ENDDO
            ENDDO ! for M
          ENDIF  ! end if (schaeffer)

          DO m = 1, smax
            CALL calc_resid_s (theta_m(1,m), a_m, b_m, m, &
               num_resid(resid_th,m), den_resid(resid_th,m), resid(resid_th,m),&
               max_resid(resid_th,m), ijk_resid(resid_th,m), zero)

            CALL under_relax_s (theta_m(1,m), a_m, b_m, m, ur_fac(8))

            CALL adjust_leq (resid(resid_th,m), leq_it(8), &
               leq_method(8), leqi, leqm)

            CALL solve_lin_eq ('Theta_m', 8, theta_m(1,m), a_m, b_m, m,&
               leqi, leqm, leq_sweep(8), leq_tol(8),  leq_pc(8), ier)
!            call out_array(Theta_m(1,m), 'Theta_m')

! Check for linear solver divergence.
            IF(ier == -2) err_l(mype) = 140

! Remove very small negative values of theta caused by leq solvers
            CALL adjust_theta (m, ier)
! large negative granular temp -> divergence
            IF (ier /= 0) err_l(mype) = 141

          ENDDO   ! end do loop (m=1,smax)
! end case(ia_2005)
! ----------------------------------------------------------------<<<


        CASE (ghd_2007)
! ---------------------------------------------------------------->>>
! do not loop over solids phases. solve for the mixture phase
          m = mmax

! initialize
          tot_no(:) = zero
          tot_sum_rs(:) = zero
          tot_eps(:) = zero

          CALL calc_nflux ()

          DO ijk = ijkstart3, ijkend3

! total number density is used for GHD theory
            cpxflux_e(ijk) = 1.5d0 * flux_ne(ijk)
            cpxflux_n(ijk) = 1.5d0 * flux_nn(ijk)
            cpxflux_t(ijk) = 1.5d0 * flux_nt(ijk)

            IF (fluid_at(ijk)) THEN
               DO l = 1,smax
                  m_pm = (pi/6.d0)*(d_p(ijk,l)**3)*ro_s(ijk,l)
                  tot_sum_rs(ijk) = tot_sum_rs(ijk) + sum_r_s(ijk,l)/m_pm
                  tot_no(ijk) = tot_no(ijk) + rop_so(ijk,l)/m_pm
                  tot_eps(ijk) = tot_eps(ijk) + ep_s(ijk,l)
               ENDDO

! calculate the source terms to be used in the a matrix and b vector
               CALL source_ghd_granular_energy (sourcelhs, &
                  sourcerhs, ijk)
               apo = 1.5d0 * tot_no(ijk)*vol(ijk)*odt
               s_p(ijk) = apo + sourcelhs + 1.5d0 *&
                  zmax(tot_sum_rs(ijk)) * vol(ijk)
               s_c(ijk) = apo*theta_mo(ijk,m) + sourcerhs + 1.50d0 *&
                  theta_m(ijk,m)*zmax(-tot_sum_rs(ijk)) * vol(ijk)
               eps(ijk) = tot_eps(ijk)
            ELSE
               eps(ijk) = zero
               s_p(ijk) = zero
               s_c(ijk) = zero
            ENDIF
          ENDDO   ! end do loop (ijk=ijkstart3,ijkend3)

! calculate the convection-diffusion terms
          CALL conv_dif_phi (theta_m(1,m), kth_s(1,m), discretize(8),&
            u_s(1,m), v_s(1,m), w_s(1,m), &
            cpxflux_e, cpxflux_n, cpxflux_t, m, a_m, b_m)

! calculate standard bc
          CALL bc_phi (theta_m(1,m), bc_theta_m(1,m), bc_thetaw_m(1,m), &
            bc_hw_theta_m(1,m), bc_c_theta_m(1,m), m, a_m, b_m)

! override bc settings if Johnson-Jackson bcs are specified
          CALL bc_theta (m, a_m, b_m)

! set the source terms in a and b matrix form
          CALL source_phi (s_p, s_c, eps, theta_m(1,m), m, a_m, b_m)

! usr source
          IF (call_usr_source(8)) CALL calc_usr_source(gran_energy,&
                                  a_m, b_m, lm=m)

          CALL calc_resid_s (theta_m(1,m), a_m, b_m, m, &
            num_resid(resid_th,m), den_resid(resid_th,m), resid(resid_th,m),&
            max_resid(resid_th,m), ijk_resid(resid_th,m), zero)

          CALL under_relax_s (theta_m(1,m), a_m, b_m, m, ur_fac(8))

          CALL adjust_leq (resid(resid_th,m), leq_it(8), leq_method(8),&
            leqi, leqm)

          CALL solve_lin_eq ('Theta_m', 8, theta_m(1,m), a_m, b_m, m, &
            leqi, leqm, leq_sweep(8), leq_tol(8),  leq_pc(8), ier)
!         call out_array(Theta_m(1,m), 'Theta_m')

! Check for linear solver divergence.
          IF(ier == -2) err_l(mype) = 140

! Remove very small negative values of theta caused by leq solvers
          CALL adjust_theta (m, ier)

! large negative granular temp -> divergence
          IF (ier /= 0) err_l(mype) = 141

! end case(ghd_2007)
! ----------------------------------------------------------------<<<

        CASE DEFAULT
! should never hit this
          WRITE (*, '(A)') 'ADJUST_THETA'
          WRITE (*, '(A,A)') 'Unknown KT_TYPE: ', kt_type
          call mfix_exit(mype)
      END SELECT   ! end selection of kt_type_enum

      call unlock_ambm

      CALL global_all_sum(err_l, err_g)
      ier = maxval(err_g)

      RETURN
      END SUBROUTINE solve_granular_energy