File: /nfs/home/0/users/jenkins/mfix.git/model/solve_energy_eq.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: SOLVE_ENERGY_EQ                                         C
4     !  Purpose: Solve energy equations                                     C
5     !                                                                      C
6     !                                                                      C
7     !  Author: M. Syamlal                                 Date: 29-APR-97  C
8     !  Reviewer:                                          Date:            C
9     !                                                                      C
10     !  Revision Number: 1                                                  C
11     !  Purpose: Eliminate energy calculations when doing DEM               C
12     !  Author: Jay Boyalakuntla                           Date: 12-Jun-04  C
13     !                                                                      C
14     !  Literature/Document References:                                     C
15     !                                                                      C
16     !  Variables referenced:                                               C
17     !  Variables modified:                                                 C
18     !  Local variables:                                                    C
19     !                                                                      C
20     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
21     
22           SUBROUTINE SOLVE_ENERGY_EQ(IER)
23     
24     !-----------------------------------------------
25     ! Modules
26     !-----------------------------------------------
27           USE param
28           USE param1
29           USE toleranc
30           USE run
31           USE physprop
32           USE geometry
33           USE fldvar
34           USE output
35           USE indices
36           USE drag
37           USE residual
38           USE ur_facs
39           USE pgcor
40           USE pscor
41           USE leqsol
42           USE bc
43           USE energy
44           USE rxns
45           Use ambm
46           Use tmp_array, S_p => ARRAY1, S_C => ARRAY2, EPs => ARRAY3
47           Use tmp_array1, VxGama => ARRAYm1
48           USE compar
49           USE discretelement
50           USE des_thermo
51           USE mflux
52           USE mpi_utility
53           USE sendrecv
54           USE ps
55           USE mms
56           USE functions
57     
58           IMPLICIT NONE
59     !-----------------------------------------------
60     ! Dummy arguments
61     !-----------------------------------------------
62     ! Error index
63           INTEGER, INTENT(INOUT) :: IER
64     !-----------------------------------------------
65     ! Local variables
66     !-----------------------------------------------
67     ! phase index
68           INTEGER :: M
69           INTEGER :: TMP_SMAX
70     !  Cp * Flux
71           DOUBLE PRECISION :: CpxFlux_E(DIMENSION_3), &
72                               CpxFlux_N(DIMENSION_3), &
73                               CpxFlux_T(DIMENSION_3)
74     ! previous time step term
75           DOUBLE PRECISION :: apo
76     ! Indices
77           INTEGER :: IJK, I, J, K
78     ! linear equation solver method and iterations
79           INTEGER :: LEQM, LEQI
80     
81     ! Arrays for storing errors:
82     ! 120 - Gas phase energy equation diverged
83     ! 121 - Solids energy equation diverged
84     ! 12x - Unclassified
85           INTEGER :: Err_l(0:numPEs-1)  ! local
86           INTEGER :: Err_g(0:numPEs-1)  ! global
87     
88     
89     ! temporary use of global arrays:
90     ! arraym1 (locally vxgama)
91     ! the volume x average gas-solids heat transfer at cell centers
92     !      DOUBLE PRECISION :: VXGAMA(DIMENSION_3, DIMENSION_M)
93     ! array1 (locally s_p)
94     ! source vector: coefficient of dependent variable
95     ! becomes part of a_m matrix; must be positive
96     !      DOUBLE PRECISION :: S_P(DIMENSION_3)
97     ! array2 (locally s_c)
98     ! source vector: constant part becomes part of b_m vector
99     !      DOUBLE PRECISION :: S_C(DIMENSION_3)
100     ! array3 (locally eps)
101     !      DOUBLE PRECISION :: EPS(DIMENSION_3)
102     
103     ! Septadiagonal matrix A_m, vector b_m
104     !      DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
105     !      DOUBLE PRECISION B_m(DIMENSION_3, 0:DIMENSION_M)
106     !-----------------------------------------------
107     
108           call lock_ambm         ! locks arrys a_m and b_m
109           call lock_tmp_array    ! locks arraym1 (locally vxgama)
110           call lock_tmp_array1   ! locks array1, array2, array3
111                                  ! (locally s_p, s_c, eps)
112     ! Initialize error flags.
113           Err_l = 0
114     
115           TMP_SMAX = SMAX
116           IF(DISCRETE_ELEMENT) THEN
117              TMP_SMAX = 0   ! Only the gas calculations are needed
118           ENDIF
119     
120     ! initializing
121           DO M = 0, TMP_SMAX
122              CALL INIT_AB_M (A_M, B_M, IJKMAX2, M, IER)
123           ENDDO
124     
125           DO IJK = IJKSTART3, IJKEND3
126              I = I_OF(IJK)
127              J = J_OF(IJK)
128              K = K_OF(IJK)
129     
130              IF (IS_ON_myPE_plus2layers(IP1(I),J,K)) THEN
131                 IF(.NOT.ADDED_MASS) THEN
132                    CpxFlux_E(IJK) = HALF * (C_pg(IJK) + C_pg(IP_OF(IJK))) * Flux_gE(IJK)
133                 ELSE
134                    CpxFlux_E(IJK) = HALF * (C_pg(IJK) + C_pg(IP_OF(IJK))) * Flux_gSE(IJK)
135                 ENDIF
136              ENDIF
137     
138              IF (IS_ON_myPE_plus2layers(I,JP1(J),K)) THEN
139                 IF(.NOT.ADDED_MASS) THEN
140                    CpxFlux_N(IJK) = HALF * (C_pg(IJK) + C_pg(JP_OF(IJK))) * Flux_gN(IJK)
141                 ELSE
142                    CpxFlux_N(IJK) = HALF * (C_pg(IJK) + C_pg(JP_OF(IJK))) * Flux_gSN(IJK)
143                 ENDIF
144              ENDIF
145     
146              IF (IS_ON_myPE_plus2layers(I,J,KP1(K))) THEN
147                 IF(.NOT.ADDED_MASS) THEN
148                    CpxFlux_T(IJK) = HALF * (C_pg(IJK) + C_pg(KP_OF(IJK))) * Flux_gT(IJK)
149                 ELSE
150                    CpxFlux_T(IJK) = HALF * (C_pg(IJK) + C_pg(KP_OF(IJK))) * Flux_gST(IJK)
151                 ENDIF
152              ENDIF
153     
154              IF (FLUID_AT(IJK)) THEN
155                 APO = ROP_GO(IJK)*C_PG(IJK)*VOL(IJK)*ODT
156                 S_P(IJK) = APO + S_RPG(IJK)*VOL(IJK)
157                 S_C(IJK) = APO*T_GO(IJK)-HOR_G(IJK)*VOL(IJK)+S_RCG(IJK)*VOL(IJK)
158                 IF(USE_MMS) S_C(IJK) = S_C(IJK) + MMS_T_G_SRC(IJK)*VOL(IJK)
159              ELSE
160                 S_P(IJK) = ZERO
161                 S_C(IJK) = ZERO
162              ENDIF
163           ENDDO
164     
165     ! Account for heat transfer between the discrete particles and the gas phase.
166           IF(DES_CONTINUUM_COUPLED) CALL DES_Hgm(S_C, S_P)
167     
168     ! calculate the convection-diffusion terms
169           CALL CONV_DIF_PHI (T_g, K_G, DISCRETIZE(6), U_G, V_G, W_G, &
170              CpxFlux_E, CpxFlux_N, CpxFlux_T, 0, A_M, B_M, IER)
171     
172     ! calculate standard bc
173           CALL BC_PHI (T_g, BC_T_G, BC_TW_G, BC_HW_T_G, BC_C_T_G, 0, A_M, B_M, IER)
174     
175     ! set the source terms in a and b matrix equation form
176           CALL SOURCE_PHI (S_P, S_C, EP_G, T_G, 0, A_M, B_M, IER)
177     
178     ! add point sources
179           IF(POINT_SOURCE) CALL POINT_SOURCE_PHI (T_g, PS_T_g, &
180              PS_CpxMFLOW_g, 0, A_M, B_M, IER)
181     
182           DO M = 1, TMP_SMAX
183              DO IJK = IJKSTART3, IJKEND3
184                 I = I_OF(IJK)
185                 J = J_OF(IJK)
186                 K = K_OF(IJK)
187     
188                 IF (IS_ON_myPE_plus2layers(IP1(I),J,K)) THEN
189                    IF(.NOT.ADDED_MASS .OR. M /= M_AM) THEN
190                       CpxFlux_E(IJK) = HALF * (C_ps(IJK,M) + C_ps(IP_OF(IJK),M)) * Flux_sE(IJK,M)
191                    ELSE   ! M=M_AM is the only phase for which virtual mass is added
192                       CpxFlux_E(IJK) = HALF * (C_ps(IJK,M) + C_ps(IP_OF(IJK),M)) * Flux_sSE(IJK)
193                    ENDIF
194                 ENDIF
195     
196                 IF (IS_ON_myPE_plus2layers(I,JP1(J),K)) THEN
197                    IF(.NOT.ADDED_MASS .OR. M /= M_AM) THEN
198                       CpxFlux_N(IJK) = HALF * (C_ps(IJK,M) + C_ps(JP_OF(IJK),M)) * Flux_sN(IJK,M)
199                    ELSE
200                       CpxFlux_N(IJK) = HALF * (C_ps(IJK,M) + C_ps(JP_OF(IJK),M)) * Flux_sSN(IJK)
201                    ENDIF
202                 ENDIF
203     
204                 IF (IS_ON_myPE_plus2layers(I,J,KP1(K))) THEN
205                    IF(.NOT.ADDED_MASS .OR. M /= M_AM) THEN
206                       CpxFlux_T(IJK) = HALF * (C_ps(IJK,M) + C_ps(KP_OF(IJK),M)) * Flux_sT(IJK,M)
207                    ELSE
208                       CpxFlux_T(IJK) = HALF * (C_ps(IJK,M) + C_ps(KP_OF(IJK),M)) * Flux_sST(IJK)
209                    ENDIF
210                 ENDIF
211     
212                 IF (FLUID_AT(IJK)) THEN
213                    APO = ROP_SO(IJK,M)*C_PS(IJK,M)*VOL(IJK)*ODT
214                    S_P(IJK) = APO + S_RPS(IJK,M)*VOL(IJK)
215                    S_C(IJK) = APO*T_SO(IJK,M) - HOR_S(IJK,M)*VOL(IJK) + &
216                       S_RCS(IJK,M)*VOL(IJK)
217                    VXGAMA(IJK,M) = GAMA_GS(IJK,M)*VOL(IJK)
218                    EPS(IJK) = EP_S(IJK,M)
219                    IF(USE_MMS) S_C(IJK) = S_C(IJK) + MMS_T_S_SRC(IJK)*VOL(IJK)
220                 ELSE
221                    S_P(IJK) = ZERO
222                    S_C(IJK) = ZERO
223                    VXGAMA(IJK,M) = ZERO
224                    EPS(IJK) = ZERO
225                    IF(USE_MMS) EPS(IJK) = EP_S(IJK,M)
226                 ENDIF
227              ENDDO   ! end do (ijk=ijkstart3,ijkend3)
228     
229     ! calculate the convection-diffusion terms
230              CALL CONV_DIF_PHI (T_s(1,M), K_S(1,M), DISCRETIZE(6), &
231                 U_S(1,M), V_S(1,M), W_S(1,M), CpxFlux_E, CpxFlux_N, &
232                 CpxFlux_T, M, A_M, B_M, IER)
233     
234     ! calculate standard bc
235              CALL BC_PHI (T_s(1,M), BC_T_S(1,M), BC_TW_S(1,M), &
236                 BC_HW_T_S(1,M), BC_C_T_S(1,M), M, A_M, B_M, IER)
237     
238     ! set the source terms in a and b matrix equation form
239              CALL SOURCE_PHI (S_P, S_C, EPS, T_S(1,M), M, A_M, B_M, IER)
240     
241     ! Add point sources.
242              IF(POINT_SOURCE) CALL POINT_SOURCE_PHI (T_s(:,M), PS_T_s(:,M),&
243                 PS_CpxMFLOW_s(:,M), M, A_M, B_M, IER)
244     
245           ENDDO   ! end do (m=1,tmp_smax)
246     
247     ! use partial elimination on interphase heat transfer term
248           IF (TMP_SMAX > 0 .AND. .NOT.USE_MMS) &
249             CALL PARTIAL_ELIM_S (T_G, T_S, VXGAMA, A_M, B_M, IER)
250     
251           CALL CALC_RESID_S (T_G, A_M, B_M, 0, NUM_RESID(RESID_T,0),&
252              DEN_RESID(RESID_T,0), RESID(RESID_T,0), MAX_RESID(RESID_T,&
253              0), IJK_RESID(RESID_T,0), ZERO, IER)
254     
255           CALL UNDER_RELAX_S (T_G, A_M, B_M, 0, UR_FAC(6), IER)
256     
257     !      call check_ab_m(a_m, b_m, 0, .false., ier)
258     !      call write_ab_m(a_m, b_m, ijkmax2, 0, ier)
259     !      write(*,*) &
260     !         resid(resid_t, 0), max_resid(resid_t, 0), &
261     !         ijk_resid(resid_t, 0)
262     
263     
264           DO M = 1, TMP_SMAX
265              CALL CALC_RESID_S (T_S(1,M), A_M, B_M, M, NUM_RESID(RESID_T,M), &
266                 DEN_RESID(RESID_T,M), RESID(RESID_T,M), MAX_RESID(&
267                 RESID_T,M), IJK_RESID(RESID_T,M), ZERO, IER)
268     
269              CALL UNDER_RELAX_S (T_S(1,M), A_M, B_M, M, UR_FAC(6), IER)
270           ENDDO
271     
272     ! set/adjust linear equation solver method and iterations
273           CALL ADJUST_LEQ(RESID(RESID_T,0), LEQ_IT(6), LEQ_METHOD(6), &
274              LEQI, LEQM, IER)
275     !      call test_lin_eq(a_m(1, -3, 0), LEQI, LEQM, LEQ_SWEEP(6), LEQ_TOL(6), LEQ_PC(6),  0, ier)
276     
277           CALL SOLVE_LIN_EQ ('T_g', 6, T_G, A_M, B_M, 0, LEQI, LEQM, &
278              LEQ_SWEEP(6), LEQ_TOL(6), LEQ_PC(6), IER)
279     ! Check for linear solver divergence.
280           IF(ier == -2) Err_l(myPE) = 120
281     
282     ! bound temperature in any fluid or flow boundary cells
283           DO IJK = IJKSTART3, IJKEND3
284              IF(.NOT.WALL_AT(IJK))&
285                 T_g(IJK) = MIN(TMAX, MAX(TMIN, T_g(IJK)))
286           ENDDO
287     
288     !      call out_array(T_g, 'T_g')
289     
290           DO M = 1, TMP_SMAX
291              CALL ADJUST_LEQ (RESID(RESID_T,M), LEQ_IT(6), LEQ_METHOD(6), &
292                 LEQI, LEQM, IER)
293     !         call test_lin_eq(a_m(1, -3, M), LEQI, LEQM, LEQ_SWEEP(6), LEQ_TOL(6), LEQ_PC(6),  0, ier)
294              CALL SOLVE_LIN_EQ ('T_s', 6, T_S(1,M), A_M, B_M, M, LEQI, &
295                 LEQM, LEQ_SWEEP(6), LEQ_TOL(6), LEQ_PC(6), IER)
296     
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, tmp_smax)
306     
307           call unlock_ambm
308           call unlock_tmp_array
309           call unlock_tmp_array1
310     
311     ! If the linear solver diverged, temperatures may take on unphysical
312     ! values. To prevent them from propogating through the domain or
313     ! causing failure in other routines, force an exit from iterate and
314     ! reduce the time step.
315           CALL global_all_sum(Err_l, Err_g)
316           IER = maxval(Err_g)
317     
318     
319           RETURN
320           END SUBROUTINE SOLVE_ENERGY_EQ
321