File: N:\mfix\model\solve_energy_eq.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: SOLVE_ENERGY_EQ                                         C
4     !  Purpose: Solve energy equations                                     C
5     !                                                                      C
6     !  Author: M. Syamlal                                 Date: 29-APR-97  C
7     !                                                                      C
8     !  Literature/Document References:                                     C
9     !                                                                      C
10     !                                                                      C
11     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
12           SUBROUTINE SOLVE_ENERGY_EQ(IER)
13     
14     ! Modules
15     !---------------------------------------------------------------------//
16           use ambm, only: a_m, b_m, lock_ambm, unlock_ambm
17           use bc, only: bc_t_g, bc_tw_g, bc_hw_t_g, bc_c_t_g
18           use bc, only: bc_t_s, bc_tw_s, bc_hw_t_s, bc_c_t_s
19           use compar, only: ijkstart3, ijkend3, mype, numpes
20           use discretelement, only: des_continuum_coupled
21           use energy, only: hor_s, s_rps, s_rcs, gama_gs
22           use energy, only: hor_g, s_rpg, s_rcg
23           use fldvar, only: ep_g, u_g, v_g, w_g, t_g, t_go, rop_go
24           use fldvar, only: ep_s, u_s, v_s, w_s, t_s, t_so, rop_so
25           use functions, only: fluid_at, wall_at
26           use functions, only: is_on_mype_plus2layers
27           use functions, only: ip_of, jp_of, kp_of
28           use geometry, only: ijkmax2, vol
29           use indices, only: i_of, j_of, k_of
30           use indices, only: ip1, jp1, kp1
31           use leqsol, only: leq_it, leq_method, leq_sweep, leq_pc, leq_tol
32           use mflux, only: flux_ge, flux_gse, flux_gn, flux_gsn
33           use mflux, only: flux_se, flux_sse, flux_sn, flux_ssn
34           use mflux, only: flux_gt, flux_gst, flux_st, flux_sst
35           use mms, only: use_mms, mms_t_g_src, mms_t_s_src
36           use mpi_utility, only: global_all_sum
37           use param, only: dimension_3, dimension_m
38           use param1, only: zero, half
39           use physprop, only: k_g, c_pg
40           use physprop, only: k_s, c_ps, smax
41           use ps, only: point_source
42           use ps, only: ps_t_g, ps_cpxmflow_g
43           use ps, only: ps_t_s, ps_cpxmflow_s
44           use residual, only: resid, max_resid, ijk_resid
45           use residual, only: num_resid, den_resid
46           use residual, only: resid_t
47           use run, only: discretize, odt, added_mass, m_am
48           use toleranc, only: tmin, tmax
49           use ur_facs, only: ur_fac
50           use usr_src, only: call_usr_source, calc_usr_source
51           use usr_src, only: gas_energy, solids_energy
52           IMPLICIT NONE
53     
54     ! Dummy arguments
55     !---------------------------------------------------------------------//
56     ! Error index
57           INTEGER, INTENT(INOUT) :: IER
58     
59     ! Local variables
60     !---------------------------------------------------------------------//
61     ! phase index
62           INTEGER :: M
63     !  Cp * Flux
64           DOUBLE PRECISION :: CpxFlux_E(DIMENSION_3), &
65                               CpxFlux_N(DIMENSION_3), &
66                               CpxFlux_T(DIMENSION_3)
67     ! previous time step term
68           DOUBLE PRECISION :: apo
69     ! Indices
70           INTEGER :: IJK, I, J, K
71     ! linear equation solver method and iterations
72           INTEGER :: LEQM, LEQI
73     
74     ! Arrays for storing errors:
75     ! 120 - Gas phase energy equation diverged
76     ! 121 - Solids energy equation diverged
77     ! 12x - Unclassified
78           INTEGER :: Err_l(0:numPEs-1)  ! local
79           INTEGER :: Err_g(0:numPEs-1)  ! global
80     
81     ! temporary use of global arrays:
82     ! arraym1 (locally vxgama)
83     ! the volume x average gas-solids heat transfer at cell centers
84     !      DOUBLE PRECISION :: VXGAMA(DIMENSION_3, DIMENSION_M)
85     ! array1 (locally s_p)
86     ! source vector: coefficient of dependent variable
87     ! becomes part of a_m matrix; must be positive
88           DOUBLE PRECISION :: S_P(DIMENSION_3)
89     ! array2 (locally s_c)
90     ! source vector: constant part becomes part of b_m vector
91           DOUBLE PRECISION :: S_C(DIMENSION_3)
92     ! array3 (locally eps)
93           DOUBLE PRECISION :: EPS(DIMENSION_3)
94     
95           DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: VXGAMA
96     
97     ! Septadiagonal matrix A_m, vector b_m
98     !      DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
99     !      DOUBLE PRECISION B_m(DIMENSION_3, 0:DIMENSION_M)
100     !---------------------------------------------------------------------//
101     
102           ALLOCATE(VXGAMA(DIMENSION_3,DIMENSION_M))
103     
104           call lock_ambm         ! locks arrys a_m and b_m
105                                  ! (locally s_p, s_c, eps)
106     ! Initialize error flags.
107           Err_l = 0
108     
109     ! initializing
110           DO M = 0, SMAX
111              CALL INIT_AB_M (A_M, B_M, IJKMAX2, M)
112           ENDDO
113     
114     ! Gas phase computations
115     ! ---------------------------------------------------------------->>>
116           DO IJK = IJKSTART3, IJKEND3
117              I = I_OF(IJK)
118              J = J_OF(IJK)
119              K = K_OF(IJK)
120     
121              IF (IS_ON_myPE_plus2layers(IP1(I),J,K)) THEN
122                 IF(.NOT.ADDED_MASS) THEN
123                    CpxFlux_E(IJK) = HALF * (C_pg(IJK) + C_pg(IP_OF(IJK))) * Flux_gE(IJK)
124                 ELSE
125                    CpxFlux_E(IJK) = HALF * (C_pg(IJK) + C_pg(IP_OF(IJK))) * Flux_gSE(IJK)
126                 ENDIF
127              ENDIF
128     
129              IF (IS_ON_myPE_plus2layers(I,JP1(J),K)) THEN
130                 IF(.NOT.ADDED_MASS) THEN
131                    CpxFlux_N(IJK) = HALF * (C_pg(IJK) + C_pg(JP_OF(IJK))) * Flux_gN(IJK)
132                 ELSE
133                    CpxFlux_N(IJK) = HALF * (C_pg(IJK) + C_pg(JP_OF(IJK))) * Flux_gSN(IJK)
134                 ENDIF
135              ENDIF
136     
137              IF (IS_ON_myPE_plus2layers(I,J,KP1(K))) THEN
138                 IF(.NOT.ADDED_MASS) THEN
139                    CpxFlux_T(IJK) = HALF * (C_pg(IJK) + C_pg(KP_OF(IJK))) * Flux_gT(IJK)
140                 ELSE
141                    CpxFlux_T(IJK) = HALF * (C_pg(IJK) + C_pg(KP_OF(IJK))) * Flux_gST(IJK)
142                 ENDIF
143              ENDIF
144     
145              IF (FLUID_AT(IJK)) THEN
146                 APO = ROP_GO(IJK)*C_PG(IJK)*VOL(IJK)*ODT
147                 S_P(IJK) = APO + S_RPG(IJK)*VOL(IJK)
148                 S_C(IJK) = APO*T_GO(IJK)-HOR_G(IJK)*VOL(IJK)+S_RCG(IJK)*VOL(IJK)
149                 IF(USE_MMS) S_C(IJK) = S_C(IJK) + MMS_T_G_SRC(IJK)*VOL(IJK)
150              ELSE
151                 S_P(IJK) = ZERO
152                 S_C(IJK) = ZERO
153              ENDIF
154           ENDDO
155     
156     ! Account for convective heat transfer between gas and DES particles.
157           IF(DES_CONTINUUM_COUPLED) CALL DES_2FLUID_CONV(S_P, S_C)
158     
159     ! calculate the convection-diffusion terms
160           CALL CONV_DIF_PHI (T_g, K_G, DISCRETIZE(6), U_G, V_G, W_G, &
161              CpxFlux_E, CpxFlux_N, CpxFlux_T, 0, A_M, B_M)
162     
163     ! calculate standard bc
164           CALL BC_PHI (T_g, BC_T_G, BC_TW_G, BC_HW_T_G, BC_C_T_G, 0, A_M, B_M)
165     
166     ! set the source terms in a and b matrix equation form
167           CALL SOURCE_PHI (S_P, S_C, EP_G, T_G, 0, A_M, B_M)
168     
169     ! add point sources
170           IF(POINT_SOURCE) CALL POINT_SOURCE_PHI (T_g, PS_T_g, &
171              PS_CpxMFLOW_g, 0, A_M, B_M)
172     ! usr source
173           IF (CALL_USR_SOURCE(6)) CALL CALC_USR_SOURCE(GAS_ENERGY, &
174                                 A_M, B_M, lM=0)
175     
176     ! Solids phases computations
177     ! ---------------------------------------------------------------->>>
178           DO M = 1, SMAX
179              DO IJK = IJKSTART3, IJKEND3
180                 I = I_OF(IJK)
181                 J = J_OF(IJK)
182                 K = K_OF(IJK)
183     
184                 IF (IS_ON_myPE_plus2layers(IP1(I),J,K)) THEN
185                    IF(.NOT.ADDED_MASS .OR. M /= M_AM) THEN
186                       CpxFlux_E(IJK) = HALF * (C_ps(IJK,M) + C_ps(IP_OF(IJK),M)) * Flux_sE(IJK,M)
187                    ELSE   ! M=M_AM is the only phase for which virtual mass is added
188                       CpxFlux_E(IJK) = HALF * (C_ps(IJK,M) + C_ps(IP_OF(IJK),M)) * Flux_sSE(IJK)
189                    ENDIF
190                 ENDIF
191     
192                 IF (IS_ON_myPE_plus2layers(I,JP1(J),K)) THEN
193                    IF(.NOT.ADDED_MASS .OR. M /= M_AM) THEN
194                       CpxFlux_N(IJK) = HALF * (C_ps(IJK,M) + C_ps(JP_OF(IJK),M)) * Flux_sN(IJK,M)
195                    ELSE
196                       CpxFlux_N(IJK) = HALF * (C_ps(IJK,M) + C_ps(JP_OF(IJK),M)) * Flux_sSN(IJK)
197                    ENDIF
198                 ENDIF
199     
200                 IF (IS_ON_myPE_plus2layers(I,J,KP1(K))) THEN
201                    IF(.NOT.ADDED_MASS .OR. M /= M_AM) THEN
202                       CpxFlux_T(IJK) = HALF * (C_ps(IJK,M) + C_ps(KP_OF(IJK),M)) * Flux_sT(IJK,M)
203                    ELSE
204                       CpxFlux_T(IJK) = HALF * (C_ps(IJK,M) + C_ps(KP_OF(IJK),M)) * Flux_sST(IJK)
205                    ENDIF
206                 ENDIF
207     
208                 IF (FLUID_AT(IJK)) THEN
209                    APO = ROP_SO(IJK,M)*C_PS(IJK,M)*VOL(IJK)*ODT
210                    S_P(IJK) = APO + S_RPS(IJK,M)*VOL(IJK)
211                    S_C(IJK) = APO*T_SO(IJK,M) - HOR_S(IJK,M)*VOL(IJK) + &
212                       S_RCS(IJK,M)*VOL(IJK)
213                    VXGAMA(IJK,M) = GAMA_GS(IJK,M)*VOL(IJK)
214                    EPS(IJK) = EP_S(IJK,M)
215                    IF(USE_MMS) S_C(IJK) = S_C(IJK) + MMS_T_S_SRC(IJK)*VOL(IJK)
216                 ELSE
217                    S_P(IJK) = ZERO
218                    S_C(IJK) = ZERO
219                    VXGAMA(IJK,M) = ZERO
220                    EPS(IJK) = ZERO
221                    IF(USE_MMS) EPS(IJK) = EP_S(IJK,M)
222                 ENDIF
223              ENDDO   ! end do (ijk=ijkstart3,ijkend3)
224     
225     ! calculate the convection-diffusion terms
226              CALL CONV_DIF_PHI (T_s(1,M), K_S(1,M), DISCRETIZE(6), &
227                 U_S(1,M), V_S(1,M), W_S(1,M), CpxFlux_E, CpxFlux_N, &
228                 CpxFlux_T, M, A_M, B_M)
229     
230     ! calculate standard bc
231              CALL BC_PHI (T_s(1,M), BC_T_S(1,M), BC_TW_S(1,M), &
232                 BC_HW_T_S(1,M), BC_C_T_S(1,M), M, A_M, B_M)
233     
234     ! set the source terms in a and b matrix equation form
235              CALL SOURCE_PHI (S_P, S_C, EPS, T_S(1,M), M, A_M, B_M)
236     
237     ! Add point sources.
238              IF(POINT_SOURCE) CALL POINT_SOURCE_PHI (T_s(:,M), PS_T_s(:,M),&
239                 PS_CpxMFLOW_s(:,M), M, A_M, B_M)
240     
241     ! usr source
242              IF (CALL_USR_SOURCE(6)) CALL CALC_USR_SOURCE(SOLIDS_ENERGY, &
243                                    A_M, B_M, lM=M)
244     
245           ENDDO   ! end do (m=1,smax)
246     
247     ! Solve gas and solids phase equations
248     ! ---------------------------------------------------------------->>>
249     ! use partial elimination on interphase heat transfer term
250           IF (SMAX > 0 .AND. .NOT.USE_MMS) &
251             CALL PARTIAL_ELIM_S (T_G, T_S, VXGAMA, A_M, B_M)
252     
253           CALL CALC_RESID_S (T_G, A_M, B_M, 0, NUM_RESID(RESID_T,0),&
254              DEN_RESID(RESID_T,0), RESID(RESID_T,0), MAX_RESID(RESID_T,&
255              0), IJK_RESID(RESID_T,0), ZERO)
256     
257           CALL UNDER_RELAX_S (T_G, A_M, B_M, 0, UR_FAC(6))
258     
259     !      call check_ab_m(a_m, b_m, 0, .false., ier)
260     !      call write_ab_m(a_m, b_m, ijkmax2, 0, ier)
261     !      write(*,*) &
262     !         resid(resid_t, 0), max_resid(resid_t, 0), &
263     !         ijk_resid(resid_t, 0)
264     
265     
266           DO M = 1, SMAX
267              CALL CALC_RESID_S (T_S(1,M), A_M, B_M, M, NUM_RESID(RESID_T,M), &
268                 DEN_RESID(RESID_T,M), RESID(RESID_T,M), MAX_RESID(&
269                 RESID_T,M), IJK_RESID(RESID_T,M), ZERO)
270     
271              CALL UNDER_RELAX_S (T_S(1,M), A_M, B_M, M, UR_FAC(6))
272           ENDDO
273     
274     ! set/adjust linear equation solver method and iterations
275           CALL ADJUST_LEQ(RESID(RESID_T,0), LEQ_IT(6), LEQ_METHOD(6), &
276              LEQI, LEQM)
277     !      call test_lin_eq(a_m(1, -3, 0), LEQI, LEQM, LEQ_SWEEP(6), LEQ_TOL(6), LEQ_PC(6),  0, ier)
278     
279           CALL SOLVE_LIN_EQ ('T_g', 6, T_G, A_M, B_M, 0, LEQI, LEQM, &
280              LEQ_SWEEP(6), LEQ_TOL(6), LEQ_PC(6), IER)
281     ! Check for linear solver divergence.
282           IF(ier == -2) Err_l(myPE) = 120
283     
284     ! bound temperature in any fluid or flow boundary cells
285           DO IJK = IJKSTART3, IJKEND3
286              IF(.NOT.WALL_AT(IJK))&
287                 T_g(IJK) = MIN(TMAX, MAX(TMIN, T_g(IJK)))
288           ENDDO
289     !      call out_array(T_g, 'T_g')
290     
291           DO M = 1, SMAX
292              CALL ADJUST_LEQ (RESID(RESID_T,M), LEQ_IT(6), LEQ_METHOD(6), &
293                 LEQI, LEQM)
294     !         call test_lin_eq(a_m(1, -3, M), LEQI, LEQM, LEQ_SWEEP(6), LEQ_TOL(6), LEQ_PC(6),  0, ier)
295              CALL SOLVE_LIN_EQ ('T_s', 6, T_S(1,M), A_M, B_M, M, LEQI, &
296                 LEQM, LEQ_SWEEP(6), LEQ_TOL(6), LEQ_PC(6), IER)
297     ! Check for linear solver divergence.
298              IF(ier == -2) Err_l(myPE) = 121
299     
300     ! bound temperature in any fluid or flow boundary cells
301              DO IJK = IJKSTART3, IJKEND3
302                 IF(.NOT.WALL_AT(IJK))&
303                    T_s(IJK, M) = MIN(TMAX, MAX(TMIN, T_s(IJK, M)))
304              ENDDO
305           ENDDO   ! end do (m=1, smax)
306     
307           call unlock_ambm
308     
309     ! If the linear solver diverged, temperatures may take on unphysical
310     ! values. To prevent them from propogating through the domain or
311     ! causing failure in other routines, force an exit from iterate and
312     ! reduce the time step.
313           CALL global_all_sum(Err_l, Err_g)
314           IER = maxval(Err_g)
315     
316           RETURN
317           END SUBROUTINE SOLVE_ENERGY_EQ
318