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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: SOLVE_GRANULAR_ENERGY                                   C
4     !  Purpose: Solve granular energy equations in matrix equation         C
5     !     form Ax=b. The center coefficient (ap) and source vector (b)     C
6     !     are negative.  The off-diagonal coefficients are positive.       C
7     !                                                                      C
8     !  Purpose: Solve granular energy equations                            C
9     !                                                                      C
10     !                                                                      C
11     !  Author: K. Agrawal                                 Date:            C
12     !  Reviewer:                                          Date:            C
13     !                                                                      C
14     !                                                                      C
15     !  Literature/Document References:                                     C
16     !                                                                      C
17     !                                                                      C
18     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
19     
20           SUBROUTINE SOLVE_GRANULAR_ENERGY(IER)
21     
22     !-----------------------------------------------
23     ! Modules
24     !-----------------------------------------------
25           USE param
26           USE param1
27           USE toleranc
28           USE run
29           USE physprop
30           USE visc_s
31           USE geometry
32           USE fldvar
33           USE constant
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 compar
48           USE mflux
49           USE mpi_utility
50           USE mms
51           USE functions
52           IMPLICIT NONE
53     !-----------------------------------------------
54     ! Dummy arguments
55     !-----------------------------------------------
56     ! Error index
57           INTEGER :: IER
58     !-----------------------------------------------
59     ! Local variables
60     !-----------------------------------------------
61     ! phase index
62           INTEGER :: M, L
63     ! Cp * Flux
64           DOUBLE PRECISION CpxFlux_E(DIMENSION_3), CpxFlux_N(DIMENSION_3), CpxFlux_T(DIMENSION_3)
65     ! previous time step term
66           DOUBLE PRECISION :: apo
67     ! source terms which appear appear in the center coefficient (lhs) and
68     ! the rhs vector
69           DOUBLE PRECISION :: sourcelhs, sourcerhs
70     ! Indices
71           INTEGER :: IJK
72     ! linear equation solver method and iterations
73           INTEGER :: LEQM, LEQI
74     ! particle mass and diameter
75           DOUBLE PRECISION :: M_PM,D_PM
76     ! total solids volume fraction
77           DOUBLE PRECISION :: TOT_EPS(DIMENSION_3)
78     ! net production of all solids phase
79           DOUBLE PRECISION :: TOT_SUM_RS(DIMENSION_3)
80     ! old value of total solids number density
81           DOUBLE PRECISION :: TOT_NO(DIMENSION_3)
82     ! small value of theta_m, 1 cm2/s2 = 1e-4 m2/s2
83           DOUBLE PRECISION :: smallTheta
84     ! temporary variable for volume x interphase transfer
85     ! coefficient between solids phases
86           DOUBLE PRECISION :: VxTC_ss(DIMENSION_3,DIMENSION_LM)
87     
88     ! Arrays for storing errors:
89     ! 140 - Linear Eq diverged
90     ! 141 - Unphysical values
91     ! 14x - Unclassified
92           INTEGER :: Err_l(0:numPEs-1)  ! local
93           INTEGER :: Err_g(0:numPEs-1)  ! global
94     
95     ! temporary use of global arrays:
96     ! array1 (locally s_p)
97     ! source vector: coefficient of dependent variable
98     ! becomes part of a_m matrix; must be positive
99     !      DOUBLE PRECISION :: S_P(DIMENSION_3)
100     ! array2 (locally s_c)
101     ! source vector: constant part becomes part of b_m vector
102     !      DOUBLE PRECISION :: S_C(DIMENSION_3)
103     ! array3 (locally eps)
104     !      DOUBLE PRECISION :: eps(DIMENSION_3)
105     
106     ! Septadiagonal matrix A_m, vector b_m
107     !      DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
108     !      DOUBLE PRECISION B_m(DIMENSION_3, 0:DIMENSION_M)
109     !-----------------------------------------------
110     
111           call lock_ambm       ! locks arrays a_m and b_m
112           call lock_tmp_array  ! locks array1,array2,array3
113                                ! (locally s_p, s_c, eps)
114     
115     ! Initialize error flags.
116           Err_l = 0
117     
118           smallTheta = (to_SI)**4 * ZERO_EP_S
119     
120           DO M = 0, MMAX
121              CALL INIT_AB_M (A_M, B_M, IJKMAX2, M, IER)
122           ENDDO
123     
124           SELECT CASE (KT_TYPE_ENUM)
125             CASE(LUN_1984, AHMADI_1995, SIMONIN_1996, GD_1999, GTSH_2012)
126     ! ---------------------------------------------------------------->>>
127               DO M = 1, SMAX
128     
129                 DO IJK = ijkstart3, ijkend3
130     
131                    IF(.NOT.ADDED_MASS .OR. M /= M_AM) THEN
132                       CpxFlux_E(IJK) = 1.5D0 * Flux_sE(IJK,M)
133                       CpxFlux_N(IJK) = 1.5D0 * Flux_sN(IJK,M)
134                       CpxFlux_T(IJK) = 1.5D0 * Flux_sT(IJK,M)
135                    ELSE
136                       CpxFlux_E(IJK) = 1.5D0 * Flux_sSE(IJK)
137                       CpxFlux_N(IJK) = 1.5D0 * Flux_sSN(IJK)
138                       CpxFlux_T(IJK) = 1.5D0 * Flux_sST(IJK)
139                    ENDIF
140     
141                    IF (FLUID_AT(IJK)) THEN
142     ! calculate the source terms to be used in the a matrix and b vector
143                       IF (KT_TYPE_ENUM == GD_1999 .OR. &
144                           KT_TYPE_ENUM == GTSH_2012) THEN
145                          CALL SOURCE_GRANULAR_ENERGY_GD(SOURCELHS, &
146                             SOURCERHS, IJK, M, IER)
147                       ELSE
148                          CALL SOURCE_GRANULAR_ENERGY(SOURCELHS, &
149                             SOURCERHS, IJK, M, IER)
150                       ENDIF
151                       APO = 1.5D0*ROP_SO(IJK,M)*VOL(IJK)*ODT
152                       S_P(IJK) = APO + SOURCELHS + 1.5D0* &
153                          ZMAX(SUM_R_S(IJK,M)) * VOL(IJK)
154                       S_C(IJK) = APO*THETA_MO(IJK,M) + SOURCERHS + 1.5d0* &
155                           THETA_M(IJK,M)*ZMAX((-SUM_R_S(IJK,M))) * VOL(IJK)
156                       EPS(IJK) = EP_S(IJK,M)
157     ! MMS Source term.
158                       IF(USE_MMS) S_C(IJK) = S_C(IJK) + &
159                          MMS_THETA_M_SRC(IJK)*VOL(IJK)
160                    ELSE
161                       EPS(IJK) = ZERO
162                       S_P(IJK) = ZERO
163                       S_C(IJK) = ZERO
164                       IF(USE_MMS) EPS(IJK) = EP_S(IJK,M)
165                    ENDIF
166                 ENDDO   ! end do loop (ijk=ijkstart3,ijkend3)
167     
168     ! calculate the convection-diffusion terms
169                 CALL CONV_DIF_PHI (THETA_M(1,M), KTH_S(1,M), DISCRETIZE(8), &
170                    U_S(1,M), V_S(1,M), W_S(1,M), &
171                    CpxFlux_E, CpxFlux_N, CpxFlux_T, M, A_M, B_M, IER)
172     
173     ! calculate standard bc
174                 CALL BC_PHI (THETA_M(1,M), BC_THETA_M(1,M), BC_THETAW_M(1,M), &
175                    BC_HW_THETA_M(1,M), BC_C_THETA_M(1,M), M, A_M, B_M, IER)
176     
177     ! override bc settings if Johnson-Jackson bcs are specified
178                 CALL BC_THETA (M, A_M, B_M, IER)
179     
180     ! set the source terms in a and b matrix form
181                 CALL SOURCE_PHI (S_P, S_C, EPS, THETA_M(1,M), M, A_M, B_M)
182     
183     ! Adjusting the values of theta_m to zero when Ep_g < EP_star
184     ! (Shaeffer, 1987). This is done here instead of calc_mu_s to
185     ! avoid convergence problems. (sof)
186                 IF (SCHAEFFER) THEN
187                    DO IJK = ijkstart3, ijkend3
188                       IF (FLUID_AT(IJK) .AND. &
189                           EP_g(IJK) .LT. EP_g_blend_start(ijk)) THEN
190     
191                          A_M(IJK,1,M) = ZERO
192                          A_M(IJK,-1,M) = ZERO
193                          A_M(IJK,2,M) = ZERO
194                          A_M(IJK,-2,M) = ZERO
195                          A_M(IJK,3,M) = ZERO
196                          A_M(IJK,-3,M) = ZERO
197                          A_M(IJK,0,M) = -ONE
198                          B_M(IJK,M) = -smallTheta
199                       ENDIF
200                    ENDDO
201                 ENDIF
202     
203                 CALL CALC_RESID_S (THETA_M(1,M), A_M, B_M, M, &
204                    NUM_RESID(RESID_TH,M), &
205                    DEN_RESID(RESID_TH,M), RESID(RESID_TH,M), &
206                    MAX_RESID(RESID_TH,M), IJK_RESID(RESID_TH,M), ZERO, IER)
207     
208                 CALL UNDER_RELAX_S (THETA_M(1,M), A_M, B_M, M, UR_FAC(8), IER)
209     
210     !            call check_ab_m(a_m, b_m, m, .true., ier)
211     !            write(*,*) resid(resid_th, m), max_resid(resid_th, m),&
212     !               ijk_resid(resid_th, m)
213     !            call write_ab_m(a_m, b_m, ijkmax2, m, ier)
214     
215     !            call test_lin_eq(ijkmax2, ijmax2, imax2, a_m(1, -3, M), &
216     !               1, DO_K, ier)
217     
218                 CALL ADJUST_LEQ (RESID(RESID_TH,M), LEQ_IT(8), LEQ_METHOD(8),&
219                    LEQI, LEQM, IER)
220     
221                 CALL SOLVE_LIN_EQ ('Theta_m', 8, THETA_M(1,M), A_M, B_M, M,&
222                    LEQI, LEQM, LEQ_SWEEP(8), LEQ_TOL(8),  LEQ_PC(8), IER)
223     !            call out_array(Theta_m(1,m), 'Theta_m')
224     
225     ! Check for linear solver divergence.
226                 IF(ier == -2) Err_l(myPE) = 140
227     
228     ! Remove very small negative values of theta caused by leq solvers
229                 CALL ADJUST_THETA (M, IER)
230     ! large negative granular temp -> divergence
231                 IF (IER /= 0) Err_l(myPE) = 141
232     
233               ENDDO   ! end do loop (m=1,smax)
234     ! end case(lun_1984, simonin_1996, ahmadi_1995, gd_1999, gtsh_2012)
235     ! ----------------------------------------------------------------<<<
236     
237             CASE(IA_2005)
238     ! ---------------------------------------------------------------->>>
239               DO M = 1, SMAX
240                 DO IJK = ijkstart3, ijkend3
241     
242     ! Skip walls where some values are undefined.
243                    IF(WALL_AT(IJK)) cycle
244     
245                    D_PM = D_P(IJK,M)
246                    M_PM = (PI/6.d0)*(D_PM**3)*RO_S(IJK,M)
247     
248     ! In Iddir & Arastoopour (2005) the granular temperature includes
249     ! mass of the particle in the definition.
250                    IF(.NOT.ADDED_MASS .OR. M /= M_AM) THEN
251                       CpxFlux_E(IJK) = (1.5D0/M_PM) * Flux_sE(IJK,M)
252                       CpxFlux_N(IJK) = (1.5D0/M_PM) * Flux_sN(IJK,M)
253                       CpxFlux_T(IJK) = (1.5D0/M_PM) * Flux_sT(IJK,M)
254                    ELSE ! in case added mass is used.
255                       CpxFlux_E(IJK) = (1.5D0/M_PM) * Flux_sSE(IJK)
256                       CpxFlux_N(IJK) = (1.5D0/M_PM) * Flux_sSN(IJK)
257                       CpxFlux_T(IJK) = (1.5D0/M_PM) * Flux_sST(IJK)
258                    ENDIF
259     
260                    IF (FLUID_AT(IJK)) THEN
261     ! calculate the source terms to be used in the a matrix and b vector
262                       CALL SOURCE_GRANULAR_ENERGY_IA(SOURCELHS, &
263                          SOURCERHS, IJK, M, IER)
264                       APO = (1.5D0/M_PM)*ROP_SO(IJK,M)*VOL(IJK)*ODT
265                       S_P(IJK) = APO + SOURCELHS + (1.5d0/M_PM)*&
266                          ZMAX(SUM_R_S(IJK,M)) * VOL(IJK)
267                       S_C(IJK) = APO*THETA_MO(IJK,M) + SOURCERHS + &
268                          (1.50d0/M_PM)*THETA_M(IJK,M)*&
269                          ZMAX((-SUM_R_S(IJK,M))) * VOL(IJK)
270                       EPS(IJK) = EP_S(IJK,M)
271                    ELSE
272                       EPS(IJK) = ZERO
273                       S_P(IJK) = ZERO
274                       S_C(IJK) = ZERO
275                    ENDIF
276                 ENDDO    ! end do loop (ijk=ijkstart3,ijkend3)
277     
278     ! calculate the convection-diffusion terms
279                 CALL CONV_DIF_PHI (THETA_M(1,M), KTH_S(1,M), DISCRETIZE(8),&
280                    U_S(1,M), V_S(1,M), W_S(1,M), &
281                    CpxFlux_E, CpxFlux_N, CpxFlux_T, M, A_M, B_M, IER)
282     
283     ! calculate standard bc
284                 CALL BC_PHI (THETA_M(1,M), BC_THETA_M(1,M), BC_THETAW_M(1,M),&
285                    BC_HW_THETA_M(1,M),BC_C_THETA_M(1,M), M, A_M, B_M, IER)
286     
287     ! override bc settings if Johnson-Jackson bcs are specified
288                 CALL BC_THETA (M, A_M, B_M, IER)
289     
290     ! set the source terms in a and b matrix form
291                 CALL SOURCE_PHI (S_P, S_C, EPS, THETA_M(1,M), M, A_M, B_M)
292               ENDDO   ! end do loop (m = 1, smax)
293     
294     ! use partial elimination on collisional dissipation term: SUM(Nip)
295     !     SUM( ED_s_ip* (Theta_p-Theta_i))
296               IF (SMAX > 1) THEN
297                 CALL CALC_VTC_SS (VXTC_SS, IER)
298                 CALL PARTIAL_ELIM_IA (THETA_M, VXTC_SS, A_M, B_M, IER)
299               ENDIF
300     
301     ! Adjusting the values of theta_m to zero when Ep_g < EP_star
302     ! (Shaeffer, 1987). This is done here instead of calc_mu_s to
303     ! avoid convergence problems. (sof)
304               IF (SCHAEFFER) THEN
305                 DO M = 1, SMAX
306                    DO IJK = ijkstart3, ijkend3
307                       IF (FLUID_AT(IJK) .AND. &
308                           EP_g(IJK) .LT. EP_g_blend_start(ijk)) THEN
309     
310                          D_PM = D_P(IJK,M)
311                          M_PM = (PI/6.d0)*(D_PM**3)*RO_S(IJK,M)
312                          A_M(IJK,1,M) = ZERO
313                          A_M(IJK,-1,M) = ZERO
314                          A_M(IJK,2,M) = ZERO
315                          A_M(IJK,-2,M) = ZERO
316                          A_M(IJK,3,M) = ZERO
317                          A_M(IJK,-3,M) = ZERO
318                          A_M(IJK,0,M) = -ONE
319     ! In Iddir & Arastoopour (2005) the granular temperature includes
320     ! mass of the particle in the definition. Systems with small
321     ! particles can give rise to smaller temperatures than the standard
322     ! zero_ep_s
323                          B_M(IJK,M) = -smallTheta*M_PM
324                       ENDIF
325                    ENDDO
326                 ENDDO ! for M
327               ENDIF  ! end if (schaeffer)
328     
329               DO M = 1, SMAX
330                 CALL CALC_RESID_S (THETA_M(1,M), A_M, B_M, M, &
331                    NUM_RESID(RESID_TH,M), DEN_RESID(RESID_TH,M), RESID(RESID_TH,M),&
332                    MAX_RESID(RESID_TH,M), IJK_RESID(RESID_TH,M), ZERO, IER)
333     
334                 CALL UNDER_RELAX_S (THETA_M(1,M), A_M, B_M, M, UR_FAC(8), IER)
335     
336                 CALL ADJUST_LEQ (RESID(RESID_TH,M), LEQ_IT(8), &
337                    LEQ_METHOD(8), LEQI, LEQM, IER)
338     
339                 CALL SOLVE_LIN_EQ ('Theta_m', 8, THETA_M(1,M), A_M, B_M, M,&
340                    LEQI, LEQM, LEQ_SWEEP(8), LEQ_TOL(8),  LEQ_PC(8), IER)
341     !            call out_array(Theta_m(1,m), 'Theta_m')
342     
343     ! Check for linear solver divergence.
344                 IF(ier == -2) Err_l(myPE) = 140
345     
346     ! Remove very small negative values of theta caused by leq solvers
347                 CALL ADJUST_THETA (M, IER)
348     ! large negative granular temp -> divergence
349                 IF (IER /= 0) Err_l(myPE) = 141
350     
351               ENDDO   ! end do loop (m=1,smax)
352     ! end case(ia_2005)
353     ! ----------------------------------------------------------------<<<
354     
355     
356             CASE (GHD_2007)
357     ! ---------------------------------------------------------------->>>
358     ! do not loop over solids phases. solve for the mixture phase
359               M = MMAX
360     
361     ! initialize
362               TOT_NO(:) = ZERO
363               TOT_SUM_RS(:) = ZERO
364               TOT_EPS(:) = ZERO
365     
366               CALL CALC_NFLUX (IER)
367     
368               DO IJK = ijkstart3, ijkend3
369     
370     ! total number density is used for GHD theory
371                 CpxFlux_E(IJK) = 1.5D0 * Flux_nE(IJK)
372                 CpxFlux_N(IJK) = 1.5D0 * Flux_nN(IJK)
373                 CpxFlux_T(IJK) = 1.5D0 * Flux_nT(IJK)
374     
375                 IF (FLUID_AT(IJK)) THEN
376                    DO L = 1,SMAX
377                       M_PM = (PI/6.d0)*(D_P(IJK,L)**3)*RO_S(IJK,L)
378                       TOT_SUM_RS(IJK) = TOT_SUM_RS(IJK) + SUM_R_S(IJK,L)/M_PM
379                       TOT_NO(IJK) = TOT_NO(IJK) + ROP_SO(IJK,L)/M_PM
380                       TOT_EPS(IJK) = TOT_EPS(IJK) + EP_S(IJK,L)
381                    ENDDO
382     
383     ! calculate the source terms to be used in the a matrix and b vector
384                    CALL SOURCE_GHD_GRANULAR_ENERGY (SOURCELHS, &
385                       SOURCERHS, IJK, IER)
386                    APO = 1.5D0 * TOT_NO(IJK)*VOL(IJK)*ODT
387                    S_P(IJK) = APO + SOURCELHS + 1.5d0 *&
388                       ZMAX(TOT_SUM_RS(IJK)) * VOL(IJK)
389                    S_C(IJK) = APO*THETA_MO(IJK,M) + SOURCERHS + 1.50d0 *&
390                       THETA_M(IJK,M)*ZMAX(-TOT_SUM_RS(IJK)) * VOL(IJK)
391                    EPS(IJK) = TOT_EPS(IJK)
392                 ELSE
393                    EPS(IJK) = ZERO
394                    S_P(IJK) = ZERO
395                    S_C(IJK) = ZERO
396                 ENDIF
397               ENDDO   ! end do loop (ijk=ijkstart3,ijkend3)
398     
399     ! calculate the convection-diffusion terms
400               CALL CONV_DIF_PHI (THETA_M(1,M), KTH_S(1,M), DISCRETIZE(8),&
401                 U_S(1,M), V_S(1,M), W_S(1,M), &
402                 CpxFlux_E, CpxFlux_N, CpxFlux_T, M, A_M, B_M, IER)
403     
404     ! calculate standard bc
405               CALL BC_PHI (THETA_M(1,M), BC_THETA_M(1,M), BC_THETAW_M(1,M), &
406                 BC_HW_THETA_M(1,M), BC_C_THETA_M(1,M), M, A_M, B_M, IER)
407     
408     ! override bc settings if Johnson-Jackson bcs are specified
409               CALL BC_THETA (M, A_M, B_M, IER)
410     
411     ! set the source terms in a and b matrix form
412               CALL SOURCE_PHI (S_P, S_C, EPS, THETA_M(1,M), M, A_M, B_M)
413     
414               CALL CALC_RESID_S (THETA_M(1,M), A_M, B_M, M, &
415                 NUM_RESID(RESID_TH,M), DEN_RESID(RESID_TH,M), RESID(RESID_TH,M),&
416                 MAX_RESID(RESID_TH,M), IJK_RESID(RESID_TH,M), ZERO, IER)
417     
418               CALL UNDER_RELAX_S (THETA_M(1,M), A_M, B_M, M, UR_FAC(8), IER)
419     
420               CALL ADJUST_LEQ (RESID(RESID_TH,M), LEQ_IT(8), LEQ_METHOD(8),&
421                 LEQI, LEQM, IER)
422     
423               CALL SOLVE_LIN_EQ ('Theta_m', 8, THETA_M(1,M), A_M, B_M, M, &
424                 LEQI, LEQM, LEQ_SWEEP(8), LEQ_TOL(8),  LEQ_PC(8), IER)
425     !         call out_array(Theta_m(1,m), 'Theta_m')
426     
427     ! Check for linear solver divergence.
428               IF(ier == -2) Err_l(myPE) = 140
429     
430     ! Remove very small negative values of theta caused by leq solvers
431               CALL ADJUST_THETA (M, IER)
432     
433     ! large negative granular temp -> divergence
434               IF (IER /= 0) Err_l(myPE) = 141
435     
436     ! end case(ghd_2007)
437     ! ----------------------------------------------------------------<<<
438     
439             CASE DEFAULT
440     ! should never hit this
441               WRITE (*, '(A)') 'ADJUST_THETA'
442               WRITE (*, '(A,A)') 'Unknown KT_TYPE: ', KT_TYPE
443               call mfix_exit(myPE)
444           END SELECT   ! end selection of kt_type_enum
445     
446           call unlock_ambm
447           call unlock_tmp_array
448     
449     
450           CALL global_all_sum(Err_l, Err_g)
451           IER = maxval(Err_g)
452     
453     
454           RETURN
455           END SUBROUTINE SOLVE_GRANULAR_ENERGY
456     
457     
458