solve_energy_eq.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: SOLVE_ENERGY_EQ                                         C
!  Purpose: Solve energy equations                                     C
!                                                                      C
!  Author: M. Syamlal                                 Date: 29-APR-97  C
!                                                                      C
!  Literature/Document References:                                     C
!                                                                      C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE solve_energy_eq(IER)

! Modules
!---------------------------------------------------------------------//
      use ambm, only: a_m, b_m, lock_ambm, unlock_ambm
      use bc, only: bc_t_g, bc_tw_g, bc_hw_t_g, bc_c_t_g
      use bc, only: bc_t_s, bc_tw_s, bc_hw_t_s, bc_c_t_s
      use compar, only: ijkstart3, ijkend3, mype, numpes
      use discretelement, only: des_continuum_coupled
      use energy, only: hor_s, s_rps, s_rcs, gama_gs
      use energy, only: hor_g, s_rpg, s_rcg
      use fldvar, only: ep_g, u_g, v_g, w_g, t_g, t_go, rop_go
      use fldvar, only: ep_s, u_s, v_s, w_s, t_s, t_so, rop_so
      use functions, only: fluid_at, wall_at
      use functions, only: is_on_mype_plus2layers
      use functions, only: ip_of, jp_of, kp_of
      use geometry, only: ijkmax2, vol
      use indices, only: i_of, j_of, k_of
      use indices, only: ip1, jp1, kp1
      use leqsol, only: leq_it, leq_method, leq_sweep, leq_pc, leq_tol
      use mflux, only: flux_ge, flux_gse, flux_gn, flux_gsn
      use mflux, only: flux_se, flux_sse, flux_sn, flux_ssn
      use mflux, only: flux_gt, flux_gst, flux_st, flux_sst
      use mms, only: use_mms, mms_t_g_src, mms_t_s_src
      use mpi_utility, only: global_all_sum
      use param, only: dimension_3, dimension_m
      use param1, only: zero, half
      use physprop, only: k_g, c_pg
      use physprop, only: k_s, c_ps, smax
      use ps, only: point_source
      use ps, only: ps_t_g, ps_cpxmflow_g
      use ps, only: ps_t_s, ps_cpxmflow_s
      use residual, only: resid, max_resid, ijk_resid
      use residual, only: num_resid, den_resid
      use residual, only: resid_t
      use run, only: discretize, odt, added_mass, m_am
      use toleranc, only: tmin, tmax
      use ur_facs, only: ur_fac
      use usr_src, only: call_usr_source, calc_usr_source
      use usr_src, only: gas_energy, solids_energy
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! Error index
      INTEGER, INTENT(INOUT) :: IER

! Local variables
!---------------------------------------------------------------------//
! phase index
      INTEGER :: M
!  Cp * Flux
      DOUBLE PRECISION :: CpxFlux_E(dimension_3), &
                          CpxFlux_N(dimension_3), &
                          CpxFlux_T(dimension_3)
! previous time step term
      DOUBLE PRECISION :: apo
! Indices
      INTEGER :: IJK, I, J, K
! linear equation solver method and iterations
      INTEGER :: LEQM, LEQI

! Arrays for storing errors:
! 120 - Gas phase energy equation diverged
! 121 - Solids energy equation diverged
! 12x - Unclassified
      INTEGER :: Err_l(0:numpes-1)  ! local
      INTEGER :: Err_g(0:numpes-1)  ! global

! temporary use of global arrays:
! arraym1 (locally vxgama)
! the volume x average gas-solids heat transfer at cell centers
!      DOUBLE PRECISION :: VXGAMA(DIMENSION_3, DIMENSION_M)
! 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)

      DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: VXGAMA

! 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)
!---------------------------------------------------------------------//

      ALLOCATE(vxgama(dimension_3,dimension_m))

      call lock_ambm         ! locks arrys a_m and b_m
                             ! (locally s_p, s_c, eps)
! Initialize error flags.
      err_l = 0

! initializing
      DO m = 0, smax
         CALL init_ab_m (a_m, b_m, ijkmax2, m)
      ENDDO

! Gas phase computations
! ---------------------------------------------------------------->>>
      DO ijk = ijkstart3, ijkend3
         i = i_of(ijk)
         j = j_of(ijk)
         k = k_of(ijk)

         IF (is_on_mype_plus2layers(ip1(i),j,k)) THEN
            IF(.NOT.added_mass) THEN
               cpxflux_e(ijk) = half * (c_pg(ijk) + c_pg(ip_of(ijk))) * flux_ge(ijk)
            ELSE
               cpxflux_e(ijk) = half * (c_pg(ijk) + c_pg(ip_of(ijk))) * flux_gse(ijk)
            ENDIF
         ENDIF

         IF (is_on_mype_plus2layers(i,jp1(j),k)) THEN
            IF(.NOT.added_mass) THEN
               cpxflux_n(ijk) = half * (c_pg(ijk) + c_pg(jp_of(ijk))) * flux_gn(ijk)
            ELSE
               cpxflux_n(ijk) = half * (c_pg(ijk) + c_pg(jp_of(ijk))) * flux_gsn(ijk)
            ENDIF
         ENDIF

         IF (is_on_mype_plus2layers(i,j,kp1(k))) THEN
            IF(.NOT.added_mass) THEN
               cpxflux_t(ijk) = half * (c_pg(ijk) + c_pg(kp_of(ijk))) * flux_gt(ijk)
            ELSE
               cpxflux_t(ijk) = half * (c_pg(ijk) + c_pg(kp_of(ijk))) * flux_gst(ijk)
            ENDIF
         ENDIF

         IF (fluid_at(ijk)) THEN
            apo = rop_go(ijk)*c_pg(ijk)*vol(ijk)*odt
            s_p(ijk) = apo + s_rpg(ijk)*vol(ijk)
            s_c(ijk) = apo*t_go(ijk)-hor_g(ijk)*vol(ijk)+s_rcg(ijk)*vol(ijk)
            IF(use_mms) s_c(ijk) = s_c(ijk) + mms_t_g_src(ijk)*vol(ijk)
         ELSE
            s_p(ijk) = zero
            s_c(ijk) = zero
         ENDIF
      ENDDO

! Account for convective heat transfer between gas and DES particles.
      IF(des_continuum_coupled) CALL des_2fluid_conv(s_p, s_c)

! calculate the convection-diffusion terms
      CALL conv_dif_phi (t_g, k_g, discretize(6), u_g, v_g, w_g, &
         cpxflux_e, cpxflux_n, cpxflux_t, 0, a_m, b_m)

! calculate standard bc
      CALL bc_phi (t_g, bc_t_g, bc_tw_g, bc_hw_t_g, bc_c_t_g, 0, a_m, b_m)

! set the source terms in a and b matrix equation form
      CALL source_phi (s_p, s_c, ep_g, t_g, 0, a_m, b_m)

! add point sources
      IF(point_source) CALL point_source_phi (t_g, ps_t_g, &
         ps_cpxmflow_g, 0, a_m, b_m)
! usr source
      IF (call_usr_source(6)) CALL calc_usr_source(gas_energy, &
                            a_m, b_m, lm=0)

! Solids phases computations
! ---------------------------------------------------------------->>>
      DO m = 1, smax
         DO ijk = ijkstart3, ijkend3
            i = i_of(ijk)
            j = j_of(ijk)
            k = k_of(ijk)

            IF (is_on_mype_plus2layers(ip1(i),j,k)) THEN
               IF(.NOT.added_mass .OR. m /= m_am) THEN
                  cpxflux_e(ijk) = half * (c_ps(ijk,m) + c_ps(ip_of(ijk),m)) * flux_se(ijk,m)
               ELSE   ! M=M_AM is the only phase for which virtual mass is added
                  cpxflux_e(ijk) = half * (c_ps(ijk,m) + c_ps(ip_of(ijk),m)) * flux_sse(ijk)
               ENDIF
            ENDIF

            IF (is_on_mype_plus2layers(i,jp1(j),k)) THEN
               IF(.NOT.added_mass .OR. m /= m_am) THEN
                  cpxflux_n(ijk) = half * (c_ps(ijk,m) + c_ps(jp_of(ijk),m)) * flux_sn(ijk,m)
               ELSE
                  cpxflux_n(ijk) = half * (c_ps(ijk,m) + c_ps(jp_of(ijk),m)) * flux_ssn(ijk)
               ENDIF
            ENDIF

            IF (is_on_mype_plus2layers(i,j,kp1(k))) THEN
               IF(.NOT.added_mass .OR. m /= m_am) THEN
                  cpxflux_t(ijk) = half * (c_ps(ijk,m) + c_ps(kp_of(ijk),m)) * flux_st(ijk,m)
               ELSE
                  cpxflux_t(ijk) = half * (c_ps(ijk,m) + c_ps(kp_of(ijk),m)) * flux_sst(ijk)
               ENDIF
            ENDIF

            IF (fluid_at(ijk)) THEN
               apo = rop_so(ijk,m)*c_ps(ijk,m)*vol(ijk)*odt
               s_p(ijk) = apo + s_rps(ijk,m)*vol(ijk)
               s_c(ijk) = apo*t_so(ijk,m) - hor_s(ijk,m)*vol(ijk) + &
                  s_rcs(ijk,m)*vol(ijk)
               vxgama(ijk,m) = gama_gs(ijk,m)*vol(ijk)
               eps(ijk) = ep_s(ijk,m)
               IF(use_mms) s_c(ijk) = s_c(ijk) + mms_t_s_src(ijk)*vol(ijk)
            ELSE
               s_p(ijk) = zero
               s_c(ijk) = zero
               vxgama(ijk,m) = zero
               eps(ijk) = zero
               IF(use_mms) eps(ijk) = ep_s(ijk,m)
            ENDIF
         ENDDO   ! end do (ijk=ijkstart3,ijkend3)

! calculate the convection-diffusion terms
         CALL conv_dif_phi (t_s(1,m), k_s(1,m), discretize(6), &
            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 (t_s(1,m), bc_t_s(1,m), bc_tw_s(1,m), &
            bc_hw_t_s(1,m), bc_c_t_s(1,m), m, a_m, b_m)

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

! Add point sources.
         IF(point_source) CALL point_source_phi (t_s(:,m), ps_t_s(:,m),&
            ps_cpxmflow_s(:,m), m, a_m, b_m)

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

      ENDDO   ! end do (m=1,smax)

! Solve gas and solids phase equations
! ---------------------------------------------------------------->>>
! use partial elimination on interphase heat transfer term
      IF (smax > 0 .AND. .NOT.use_mms) &
        CALL partial_elim_s (t_g, t_s, vxgama, a_m, b_m)

      CALL calc_resid_s (t_g, a_m, b_m, 0, num_resid(resid_t,0),&
         den_resid(resid_t,0), resid(resid_t,0), max_resid(resid_t,&
         0), ijk_resid(resid_t,0), zero)

      CALL under_relax_s (t_g, a_m, b_m, 0, ur_fac(6))

!      call check_ab_m(a_m, b_m, 0, .false., ier)
!      call write_ab_m(a_m, b_m, ijkmax2, 0, ier)
!      write(*,*) &
!         resid(resid_t, 0), max_resid(resid_t, 0), &
!         ijk_resid(resid_t, 0)


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

         CALL under_relax_s (t_s(1,m), a_m, b_m, m, ur_fac(6))
      ENDDO

! set/adjust linear equation solver method and iterations
      CALL adjust_leq(resid(resid_t,0), leq_it(6), leq_method(6), &
         leqi, leqm)
!      call test_lin_eq(a_m(1, -3, 0), LEQI, LEQM, LEQ_SWEEP(6), LEQ_TOL(6), LEQ_PC(6),  0, ier)

      CALL solve_lin_eq ('T_g', 6, t_g, a_m, b_m, 0, leqi, leqm, &
         leq_sweep(6), leq_tol(6), leq_pc(6), ier)
! Check for linear solver divergence.
      IF(ier == -2) err_l(mype) = 120

! bound temperature in any fluid or flow boundary cells
      DO ijk = ijkstart3, ijkend3
         IF(.NOT.wall_at(ijk))&
            t_g(ijk) = min(tmax, max(tmin, t_g(ijk)))
      ENDDO
!      call out_array(T_g, 'T_g')

      DO m = 1, smax
         CALL adjust_leq (resid(resid_t,m), leq_it(6), leq_method(6), &
            leqi, leqm)
!         call test_lin_eq(a_m(1, -3, M), LEQI, LEQM, LEQ_SWEEP(6), LEQ_TOL(6), LEQ_PC(6),  0, ier)
         CALL solve_lin_eq ('T_s', 6, t_s(1,m), a_m, b_m, m, leqi, &
            leqm, leq_sweep(6), leq_tol(6), leq_pc(6), ier)
! Check for linear solver divergence.
         IF(ier == -2) err_l(mype) = 121

! bound temperature in any fluid or flow boundary cells
         DO ijk = ijkstart3, ijkend3
            IF(.NOT.wall_at(ijk))&
               t_s(ijk, m) = min(tmax, max(tmin, t_s(ijk, m)))
         ENDDO
      ENDDO   ! end do (m=1, smax)

      call unlock_ambm

! If the linear solver diverged, temperatures may take on unphysical
! values. To prevent them from propogating through the domain or
! causing failure in other routines, force an exit from iterate and
! reduce the time step.
      CALL global_all_sum(err_l, err_g)
      ier = maxval(err_g)

      RETURN
      END SUBROUTINE solve_energy_eq