File: N:\mfix\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     !  Author: K. Agrawal                                 Date:            C
9     !                                                                      C
10     !  Literature/Document References:                                     C
11     !                                                                      C
12     !                                                                      C
13     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
14           SUBROUTINE SOLVE_GRANULAR_ENERGY(IER)
15     
16     ! Modules
17     !---------------------------------------------------------------------//
18           use ambm, only: a_m, b_m, lock_ambm, unlock_ambm
19           use bc, only: bc_theta_m, bc_thetaw_m, bc_hw_theta_m, bc_c_theta_m
20           use compar, only: ijkstart3, ijkend3, mype, numpes
21           use constant, only: to_si, pi
22           use exit, only: mfix_exit
23           use fldvar, only: ep_g, theta_m, theta_mo
24           use fldvar, only: ep_s, u_s, v_s, w_s, rop_so
25           use fldvar, only: ro_s, d_p
26           use functions, only: fluid_at, wall_at, zmax
27           use geometry, only: ijkmax2, vol
28           use leqsol, only: leq_it, leq_method, leq_sweep, leq_pc, leq_tol
29           use mflux, only: flux_se, flux_sn, flux_st
30           use mflux, only: flux_sse, flux_ssn, flux_sst
31           use mflux, only: flux_ne, flux_nn, flux_nt
32           use mms, only: use_mms, mms_theta_m_src
33           use mpi_utility, only: global_all_sum
34           use param, only: dimension_3, dimension_m
35           use param1, only: zero, one, dimension_lm
36           use physprop, only: kth_s
37           use physprop, only: smax, mmax
38           use residual, only: resid, max_resid, ijk_resid
39           use residual, only: num_resid, den_resid
40           use residual, only: resid_th
41           use run, only: discretize, odt, added_mass, m_am, schaeffer
42           use run, only: kt_type_enum, lun_1984, ahmadi_1995, simonin_1996
43           use run, only: kt_type, gd_1999, gtsh_2012, ghd_2007, ia_2005
44           use rxns, only: sum_r_s
45           use toleranc, only: zero_ep_s 
46           use ur_facs, only: ur_fac
47           use usr_src, only: call_usr_source, calc_usr_source
48           use usr_src, only: gran_energy
49           use visc_s, only: ep_g_blend_start
50           IMPLICIT NONE
51     
52     ! Dummy arguments
53     !---------------------------------------------------------------------//
54     ! Error index
55           INTEGER :: IER
56     
57     ! Local variables
58     !---------------------------------------------------------------------//
59     ! phase index
60           INTEGER :: M, L
61     ! Cp * Flux
62           DOUBLE PRECISION CpxFlux_E(DIMENSION_3), CpxFlux_N(DIMENSION_3), CpxFlux_T(DIMENSION_3)
63     ! previous time step term
64           DOUBLE PRECISION :: apo
65     ! source terms which appear appear in the center coefficient (lhs) and
66     ! the rhs vector
67           DOUBLE PRECISION :: sourcelhs, sourcerhs
68     ! Indices
69           INTEGER :: IJK
70     ! linear equation solver method and iterations
71           INTEGER :: LEQM, LEQI
72     ! particle mass and diameter
73           DOUBLE PRECISION :: M_PM,D_PM
74     ! total solids volume fraction
75           DOUBLE PRECISION :: TOT_EPS(DIMENSION_3)
76     ! net production of all solids phase
77           DOUBLE PRECISION :: TOT_SUM_RS(DIMENSION_3)
78     ! old value of total solids number density
79           DOUBLE PRECISION :: TOT_NO(DIMENSION_3)
80     ! small value of theta_m, 1 cm2/s2 = 1e-4 m2/s2
81           DOUBLE PRECISION :: smallTheta
82     ! temporary variable for volume x interphase transfer
83     ! coefficient between solids phases
84           DOUBLE PRECISION :: VxTC_ss(DIMENSION_3,DIMENSION_LM)
85     
86     ! Arrays for storing errors:
87     ! 140 - Linear Eq diverged
88     ! 141 - Unphysical values
89     ! 14x - Unclassified
90           INTEGER :: Err_l(0:numPEs-1)  ! local
91           INTEGER :: Err_g(0:numPEs-1)  ! global
92     
93     ! temporary use of global arrays:
94     ! array1 (locally s_p)
95     ! source vector: coefficient of dependent variable
96     ! becomes part of a_m matrix; must be positive
97           DOUBLE PRECISION :: S_P(DIMENSION_3)
98     ! array2 (locally s_c)
99     ! source vector: constant part becomes part of b_m vector
100           DOUBLE PRECISION :: S_C(DIMENSION_3)
101     ! array3 (locally eps)
102           DOUBLE PRECISION :: eps(DIMENSION_3)
103     
104     ! Septadiagonal matrix A_m, vector b_m
105     !      DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
106     !      DOUBLE PRECISION B_m(DIMENSION_3, 0:DIMENSION_M)
107     !---------------------------------------------------------------------//
108     
109           call lock_ambm       ! locks arrays a_m and b_m
110     
111     ! Initialize error flags.
112           Err_l = 0
113     
114           smallTheta = (to_SI)**4 * ZERO_EP_S
115     
116           DO M = 0, MMAX
117              CALL INIT_AB_M (A_M, B_M, IJKMAX2, M)
118           ENDDO
119     
120           SELECT CASE (KT_TYPE_ENUM)
121             CASE(LUN_1984, AHMADI_1995, SIMONIN_1996, GD_1999, GTSH_2012)
122     ! ---------------------------------------------------------------->>>
123               DO M = 1, SMAX
124     
125                 DO IJK = ijkstart3, ijkend3
126     
127                    IF(.NOT.ADDED_MASS .OR. M /= M_AM) THEN
128                       CpxFlux_E(IJK) = 1.5D0 * Flux_sE(IJK,M)
129                       CpxFlux_N(IJK) = 1.5D0 * Flux_sN(IJK,M)
130                       CpxFlux_T(IJK) = 1.5D0 * Flux_sT(IJK,M)
131                    ELSE
132                       CpxFlux_E(IJK) = 1.5D0 * Flux_sSE(IJK)
133                       CpxFlux_N(IJK) = 1.5D0 * Flux_sSN(IJK)
134                       CpxFlux_T(IJK) = 1.5D0 * Flux_sST(IJK)
135                    ENDIF
136     
137                    IF (FLUID_AT(IJK)) THEN
138     ! calculate the source terms to be used in the a matrix and b vector
139                       IF (KT_TYPE_ENUM == GD_1999 .OR. &
140                           KT_TYPE_ENUM == GTSH_2012) THEN
141                          CALL SOURCE_GRANULAR_ENERGY_GD(SOURCELHS, &
142                             SOURCERHS, IJK, M)
143                       ELSE
144                          CALL SOURCE_GRANULAR_ENERGY(SOURCELHS, &
145                             SOURCERHS, IJK, M)
146                       ENDIF
147                       APO = 1.5D0*ROP_SO(IJK,M)*VOL(IJK)*ODT
148                       S_P(IJK) = APO + SOURCELHS + 1.5D0* &
149                          ZMAX(SUM_R_S(IJK,M)) * VOL(IJK)
150                       S_C(IJK) = APO*THETA_MO(IJK,M) + SOURCERHS + 1.5d0* &
151                           THETA_M(IJK,M)*ZMAX((-SUM_R_S(IJK,M))) * VOL(IJK)
152                       EPS(IJK) = EP_S(IJK,M)
153     ! MMS Source term.
154                       IF(USE_MMS) S_C(IJK) = S_C(IJK) + &
155                          MMS_THETA_M_SRC(IJK)*VOL(IJK)
156                    ELSE
157                       EPS(IJK) = ZERO
158                       S_P(IJK) = ZERO
159                       S_C(IJK) = ZERO
160                       IF(USE_MMS) EPS(IJK) = EP_S(IJK,M)
161                    ENDIF
162                 ENDDO   ! end do loop (ijk=ijkstart3,ijkend3)
163     
164     ! calculate the convection-diffusion terms
165                 CALL CONV_DIF_PHI (THETA_M(1,M), KTH_S(1,M), DISCRETIZE(8), &
166                    U_S(1,M), V_S(1,M), W_S(1,M), &
167                    CpxFlux_E, CpxFlux_N, CpxFlux_T, M, A_M, B_M)
168     
169     ! calculate standard bc
170                 CALL BC_PHI (THETA_M(1,M), BC_THETA_M(1,M), BC_THETAW_M(1,M), &
171                    BC_HW_THETA_M(1,M), BC_C_THETA_M(1,M), M, A_M, B_M)
172     
173     ! override bc settings if Johnson-Jackson bcs are specified
174                 CALL BC_THETA (M, A_M, B_M)
175     
176     ! set the source terms in a and b matrix form
177                 CALL SOURCE_PHI (S_P, S_C, EPS, THETA_M(1,M), M, A_M, B_M)
178     
179     ! usr source
180                 IF (CALL_USR_SOURCE(8)) CALL CALC_USR_SOURCE(GRAN_ENERGY,&
181                                       A_M, B_M, lM=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)
207     
208                 CALL UNDER_RELAX_S (THETA_M(1,M), A_M, B_M, M, UR_FAC(8))
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)
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)
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)
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)
286     
287     ! override bc settings if Johnson-Jackson bcs are specified
288                 CALL BC_THETA (M, A_M, B_M)
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     
293     ! usr source
294                 IF (CALL_USR_SOURCE(8)) CALL CALC_USR_SOURCE(GRAN_ENERGY,&
295                                       A_M, B_M, lM=M)
296               ENDDO   ! end do loop (m = 1, smax)
297     
298     ! use partial elimination on collisional dissipation term: SUM(Nip)
299     !     SUM( ED_s_ip* (Theta_p-Theta_i))
300               IF (SMAX > 1) THEN
301                 CALL CALC_VTC_SS (VXTC_SS)
302                 CALL PARTIAL_ELIM_IA (THETA_M, VXTC_SS, A_M, B_M)
303               ENDIF
304     
305     ! Adjusting the values of theta_m to zero when Ep_g < EP_star
306     ! (Shaeffer, 1987). This is done here instead of calc_mu_s to
307     ! avoid convergence problems. (sof)
308               IF (SCHAEFFER) THEN
309                 DO M = 1, SMAX
310                    DO IJK = ijkstart3, ijkend3
311                       IF (FLUID_AT(IJK) .AND. &
312                           EP_g(IJK) .LT. EP_g_blend_start(ijk)) THEN
313     
314                          D_PM = D_P(IJK,M)
315                          M_PM = (PI/6.d0)*(D_PM**3)*RO_S(IJK,M)
316                          A_M(IJK,1,M) = ZERO
317                          A_M(IJK,-1,M) = ZERO
318                          A_M(IJK,2,M) = ZERO
319                          A_M(IJK,-2,M) = ZERO
320                          A_M(IJK,3,M) = ZERO
321                          A_M(IJK,-3,M) = ZERO
322                          A_M(IJK,0,M) = -ONE
323     ! In Iddir & Arastoopour (2005) the granular temperature includes
324     ! mass of the particle in the definition. Systems with small
325     ! particles can give rise to smaller temperatures than the standard
326     ! zero_ep_s
327                          B_M(IJK,M) = -smallTheta*M_PM
328                       ENDIF
329                    ENDDO
330                 ENDDO ! for M
331               ENDIF  ! end if (schaeffer)
332     
333               DO M = 1, SMAX
334                 CALL CALC_RESID_S (THETA_M(1,M), A_M, B_M, M, &
335                    NUM_RESID(RESID_TH,M), DEN_RESID(RESID_TH,M), RESID(RESID_TH,M),&
336                    MAX_RESID(RESID_TH,M), IJK_RESID(RESID_TH,M), ZERO)
337     
338                 CALL UNDER_RELAX_S (THETA_M(1,M), A_M, B_M, M, UR_FAC(8))
339     
340                 CALL ADJUST_LEQ (RESID(RESID_TH,M), LEQ_IT(8), &
341                    LEQ_METHOD(8), LEQI, LEQM)
342     
343                 CALL SOLVE_LIN_EQ ('Theta_m', 8, THETA_M(1,M), A_M, B_M, M,&
344                    LEQI, LEQM, LEQ_SWEEP(8), LEQ_TOL(8),  LEQ_PC(8), IER)
345     !            call out_array(Theta_m(1,m), 'Theta_m')
346     
347     ! Check for linear solver divergence.
348                 IF(ier == -2) Err_l(myPE) = 140
349     
350     ! Remove very small negative values of theta caused by leq solvers
351                 CALL ADJUST_THETA (M, IER)
352     ! large negative granular temp -> divergence
353                 IF (IER /= 0) Err_l(myPE) = 141
354     
355               ENDDO   ! end do loop (m=1,smax)
356     ! end case(ia_2005)
357     ! ----------------------------------------------------------------<<<
358     
359     
360             CASE (GHD_2007)
361     ! ---------------------------------------------------------------->>>
362     ! do not loop over solids phases. solve for the mixture phase
363               M = MMAX
364     
365     ! initialize
366               TOT_NO(:) = ZERO
367               TOT_SUM_RS(:) = ZERO
368               TOT_EPS(:) = ZERO
369     
370               CALL CALC_NFLUX ()
371     
372               DO IJK = ijkstart3, ijkend3
373     
374     ! total number density is used for GHD theory
375                 CpxFlux_E(IJK) = 1.5D0 * Flux_nE(IJK)
376                 CpxFlux_N(IJK) = 1.5D0 * Flux_nN(IJK)
377                 CpxFlux_T(IJK) = 1.5D0 * Flux_nT(IJK)
378     
379                 IF (FLUID_AT(IJK)) THEN
380                    DO L = 1,SMAX
381                       M_PM = (PI/6.d0)*(D_P(IJK,L)**3)*RO_S(IJK,L)
382                       TOT_SUM_RS(IJK) = TOT_SUM_RS(IJK) + SUM_R_S(IJK,L)/M_PM
383                       TOT_NO(IJK) = TOT_NO(IJK) + ROP_SO(IJK,L)/M_PM
384                       TOT_EPS(IJK) = TOT_EPS(IJK) + EP_S(IJK,L)
385                    ENDDO
386     
387     ! calculate the source terms to be used in the a matrix and b vector
388                    CALL SOURCE_GHD_GRANULAR_ENERGY (SOURCELHS, &
389                       SOURCERHS, IJK)
390                    APO = 1.5D0 * TOT_NO(IJK)*VOL(IJK)*ODT
391                    S_P(IJK) = APO + SOURCELHS + 1.5d0 *&
392                       ZMAX(TOT_SUM_RS(IJK)) * VOL(IJK)
393                    S_C(IJK) = APO*THETA_MO(IJK,M) + SOURCERHS + 1.50d0 *&
394                       THETA_M(IJK,M)*ZMAX(-TOT_SUM_RS(IJK)) * VOL(IJK)
395                    EPS(IJK) = TOT_EPS(IJK)
396                 ELSE
397                    EPS(IJK) = ZERO
398                    S_P(IJK) = ZERO
399                    S_C(IJK) = ZERO
400                 ENDIF
401               ENDDO   ! end do loop (ijk=ijkstart3,ijkend3)
402     
403     ! calculate the convection-diffusion terms
404               CALL CONV_DIF_PHI (THETA_M(1,M), KTH_S(1,M), DISCRETIZE(8),&
405                 U_S(1,M), V_S(1,M), W_S(1,M), &
406                 CpxFlux_E, CpxFlux_N, CpxFlux_T, M, A_M, B_M)
407     
408     ! calculate standard bc
409               CALL BC_PHI (THETA_M(1,M), BC_THETA_M(1,M), BC_THETAW_M(1,M), &
410                 BC_HW_THETA_M(1,M), BC_C_THETA_M(1,M), M, A_M, B_M)
411     
412     ! override bc settings if Johnson-Jackson bcs are specified
413               CALL BC_THETA (M, A_M, B_M)
414     
415     ! set the source terms in a and b matrix form
416               CALL SOURCE_PHI (S_P, S_C, EPS, THETA_M(1,M), M, A_M, B_M)
417     
418     ! usr source
419               IF (CALL_USR_SOURCE(8)) CALL CALC_USR_SOURCE(GRAN_ENERGY,&
420                                       A_M, B_M, lM=M)
421     
422               CALL CALC_RESID_S (THETA_M(1,M), A_M, B_M, M, &
423                 NUM_RESID(RESID_TH,M), DEN_RESID(RESID_TH,M), RESID(RESID_TH,M),&
424                 MAX_RESID(RESID_TH,M), IJK_RESID(RESID_TH,M), ZERO)
425     
426               CALL UNDER_RELAX_S (THETA_M(1,M), A_M, B_M, M, UR_FAC(8))
427     
428               CALL ADJUST_LEQ (RESID(RESID_TH,M), LEQ_IT(8), LEQ_METHOD(8),&
429                 LEQI, LEQM)
430     
431               CALL SOLVE_LIN_EQ ('Theta_m', 8, THETA_M(1,M), A_M, B_M, M, &
432                 LEQI, LEQM, LEQ_SWEEP(8), LEQ_TOL(8),  LEQ_PC(8), IER)
433     !         call out_array(Theta_m(1,m), 'Theta_m')
434     
435     ! Check for linear solver divergence.
436               IF(ier == -2) Err_l(myPE) = 140
437     
438     ! Remove very small negative values of theta caused by leq solvers
439               CALL ADJUST_THETA (M, IER)
440     
441     ! large negative granular temp -> divergence
442               IF (IER /= 0) Err_l(myPE) = 141
443     
444     ! end case(ghd_2007)
445     ! ----------------------------------------------------------------<<<
446     
447             CASE DEFAULT
448     ! should never hit this
449               WRITE (*, '(A)') 'ADJUST_THETA'
450               WRITE (*, '(A,A)') 'Unknown KT_TYPE: ', KT_TYPE
451               call mfix_exit(myPE)
452           END SELECT   ! end selection of kt_type_enum
453     
454           call unlock_ambm
455     
456           CALL global_all_sum(Err_l, Err_g)
457           IER = maxval(Err_g)
458     
459           RETURN
460           END SUBROUTINE SOLVE_GRANULAR_ENERGY
461