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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: SOLVE_SPECIES_EQ                                        C
4     !  Purpose: Solve species mass balance 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: M. Syamlal                                 Date: 11-FEB-98  C
9     !  Reviewer:                                          Date:            C
10     !                                                                      C
11     !  Literature/Document References:                                     C
12     !                                                                      C
13     !  Variables referenced:                                               C
14     !  Variables modified:                                                 C
15     !  Local variables:                                                    C
16     !                                                                      C
17     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
18     
19           SUBROUTINE SOLVE_SPECIES_EQ(IER)
20     
21     !-----------------------------------------------
22     ! Modules
23     !-----------------------------------------------
24           USE ChiScheme
25           USE bc
26           USE compar
27           USE drag
28           USE energy
29           USE fldvar
30           USE geometry
31           USE indices
32           USE leqsol
33           USE matrix
34           USE mflux
35           USE mpi_utility
36           USE output
37           USE param
38           USE param1
39           USE pgcor
40           USE physprop
41           USE pscor
42           USE residual
43           USE run
44           USE rxns
45           USE sendrecv
46           USE toleranc
47           USE ur_facs
48           Use ambm
49           Use tmp_array, S_p => Array1, S_c => Array2, EPs => Array3, VxGama => Array4
50           use functions
51           use ps
52           use toleranc
53     
54           IMPLICIT NONE
55     !-----------------------------------------------
56     ! Dummy arguments
57     !-----------------------------------------------
58     ! Error index
59           INTEGER, INTENT(INOUT) :: IER
60     !-----------------------------------------------
61     ! Local variables
62     !-----------------------------------------------
63     ! phase index
64           INTEGER :: M
65     ! species index
66           INTEGER :: LN
67     ! previous time step term
68           DOUBLE PRECISION :: apo
69     ! Indices
70           INTEGER :: IJK
71     ! linear equation solver method and iterations
72           INTEGER :: LEQM, LEQI
73     ! tmp array to pass to set_chi
74           DOUBLE PRECISION :: X_s_temp(DIMENSION_3, DIMENSION_N_s)
75     
76     ! Arrays for storing errors:
77     ! 130 - Gas phase species equation diverged
78     ! 131 - Solids phase species equation diverged
79     ! 13x - Unclassified
80           INTEGER :: Err_l(0:numPEs-1)  ! local
81           INTEGER :: Err_g(0:numPEs-1)  ! global
82     
83     ! temporary use of global arrays:
84     ! array1 (locally s_p)
85     ! source lhs: coefficient of dependent variable
86     ! becomes part of a_m matrix; must be positive
87     !      DOUBLE PRECISION :: S_P(DIMENSION_3)
88     ! array2 (locally s_c)
89     ! source rhs vector: constant part becomes part of b_m vector
90     !      DOUBLE PRECISION :: S_C(DIMENSION_3)
91     ! array3 (locally eps)
92     ! alias for solids volume fraction
93     !      DOUBLE PRECISION :: eps(DIMENSION_3)
94     ! array4 (locally vxgama)
95     !      DOUBLE PRECISION :: vxgama(DIMENSION_3)
96     ! Septadiagonal matrix A_m, vector b_m
97     !      DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
98     !      DOUBLE PRECISION B_m(DIMENSION_3, 0:DIMENSION_M)
99     
100     !-----------------------------------------------
101     ! External functions
102     !-----------------------------------------------
103           DOUBLE PRECISION , EXTERNAL :: Check_conservation
104     !-----------------------------------------------
105     
106           call lock_ambm       ! locks arrys a_m and b_m
107           call lock_tmp_array  ! locks array1,array2,array3,array4
108                                ! (locally s_p, s_c, eps, vxgama)
109     
110     ! Initialize error flag.
111           Err_l = 0
112     
113     ! Fluid phase species mass balance equations
114     ! ---------------------------------------------------------------->>>
115           IF (SPECIES_EQ(0)) THEN
116              IF(chi_scheme) call set_chi(DISCRETIZE(7), X_g, NMAX(0), &
117                                          U_g, V_g, W_g)
118     
119     ! looping over species
120              DO LN = 1, NMAX(0)
121                 CALL INIT_AB_M (A_M, B_M, IJKMAX2, 0, IER)
122     !!$omp    parallel do private(IJK, APO)
123                 DO IJK = ijkstart3, ijkend3
124                    IF (FLUID_AT(IJK)) THEN
125     ! calculate the source terms to be used in the a matrix and b vector
126                        APO = ROP_GO(IJK)*VOL(IJK)*ODT
127                        S_P(IJK) = APO+&
128                           (ZMAX(SUM_R_G(IJK))+ROX_GC(IJK,LN))*VOL(IJK)
129                        S_C(IJK) = APO*X_GO(IJK,LN) + &
130                           X_G(IJK,LN)*ZMAX((-SUM_R_G(IJK)))*VOL(IJK) +&
131                           R_GP(IJK,LN)*VOL(IJK)
132                    ELSE
133                       S_P(IJK) = ZERO
134                       S_C(IJK) = ZERO
135                    ENDIF
136                 ENDDO
137     
138     ! calculate the convection-diffusion terms
139                 IF(.NOT.ADDED_MASS) THEN
140                    CALL CONV_DIF_PHI (X_G(1,LN), DIF_G(1,LN), &
141                       DISCRETIZE(7), U_G, V_G, W_G, &
142                       Flux_gE, Flux_gN, Flux_gT, 0, A_M, B_M, IER)
143                 ELSE
144                    CALL CONV_DIF_PHI (X_G(1,LN), DIF_G(1,LN),&
145                       DISCRETIZE(7), U_G, V_G, W_G, &
146                       Flux_gSE, Flux_gSN, Flux_gST, 0, A_M, B_M, IER)
147                 ENDIF
148     
149     ! calculate standard bc
150                 CALL BC_PHI (X_G(1,LN), BC_X_G(1,LN), BC_XW_G(1,LN), &
151                    BC_HW_X_G(1,LN), BC_C_X_G(1,LN), 0, A_M, B_M, IER)
152     
153     ! set the source terms in a and b matrix form
154                 CALL SOURCE_PHI (S_P, S_C, EP_G, X_G(1,LN), 0, A_M, B_M)
155     
156     ! Add point sources.
157                 IF(POINT_SOURCE) CALL POINT_SOURCE_PHI (X_G(1,LN), &
158                    PS_X_G(:,LN), PS_MASSFLOW_G, 0, A_M, B_M, IER)
159     
160                 CALL CALC_RESID_S (X_G(1,LN), A_M, B_M, 0, &
161                    NUM_RESID(RESID_X+(LN-1),0), &
162                    DEN_RESID(RESID_X+(LN-1),0), RESID(RESID_X+(LN-1),0), &
163                    MAX_RESID(RESID_X+(LN-1),0), IJK_RESID(RESID_X+(LN-1),0), &
164                    ZERO_X_GS, IER)
165     
166                 CALL UNDER_RELAX_S (X_G(1,LN), A_M, B_M, 0, UR_FAC(7), IER)
167     
168                 CALL ADJUST_LEQ (RESID(RESID_X+(LN-1),0), LEQ_IT(7), &
169                    LEQ_METHOD(7), LEQI, LEQM, IER)
170     
171     ! solve the phi equation
172                 CALL SOLVE_LIN_EQ ('X_g', 7, X_G(1,LN), A_M, B_M, 0, &
173                    LEQI, LEQM, LEQ_SWEEP(7), LEQ_TOL(7), LEQ_PC(7), IER)
174     
175     ! Check for linear solver divergence.
176                    IF(ier == -2) Err_l(myPE) = 130
177     
178                 CALL BOUND_X (X_G(1,LN), IJKMAX2, IER)
179     
180              ENDDO    ! end do loop (ln = 1, nmax(0)
181              IF(chi_scheme) call unset_chi()
182           ENDIF
183     ! end fluid phase species equations
184     ! ----------------------------------------------------------------<<<
185     
186     ! Solids phase species balance equations
187     ! ---------------------------------------------------------------->>>
188           DO M = 1, SMAX
189              IF (SPECIES_EQ(M)) THEN
190                 IF(chi_scheme) THEN
191                    DO LN = 1, NMAX(M)
192                       DO IJK = ijkstart3, ijkend3
193                          X_S_temp(IJK, LN) = X_S(IJK,M,LN)
194                       ENDDO
195                     ENDDO
196                   call set_chi(DISCRETIZE(7), X_S_temp, NMAX(M), &
197                                U_S(1,M), V_S(1,M), W_S(1,M))
198                 ENDIF ! for chi_scheme
199     
200                 DO LN = 1, NMAX(M)
201                    CALL INIT_AB_M (A_M, B_M, IJKMAX2, M, IER)
202     
203     !!$omp    parallel do private(IJK, APO)
204                    DO IJK = ijkstart3, ijkend3
205                       IF (FLUID_AT(IJK)) THEN
206                         APO = ROP_SO(IJK,M)*VOL(IJK)*ODT
207                         S_P(IJK) = APO + &
208                            (ZMAX(SUM_R_S(IJK,M))+ROX_SC(IJK,M,LN))*VOL(IJK)
209                         S_C(IJK) = APO*X_SO(IJK,M,LN) + &
210                            X_S(IJK,M,LN)*ZMAX((-SUM_R_S(IJK,M)))*VOL(IJK) + &
211                            R_SP(IJK,M,LN)*VOL(IJK)
212                         EPS(IJK) = EP_S(IJK,M)
213                       ELSE
214                          S_P(IJK) = ZERO
215                          S_C(IJK) = ZERO
216                          EPS(IJK) = ZERO
217                       ENDIF
218                    ENDDO
219     
220                    IF(.NOT.ADDED_MASS .OR. M /= M_AM) THEN
221                       CALL CONV_DIF_PHI (X_S(1,M,LN), DIF_S(1,M,LN), &
222                         DISCRETIZE(7), U_S(1,M), V_S(1,M), W_S(1,M), &
223                         Flux_sE(1,M), Flux_sN(1,M), Flux_sT(1,M), M, A_M, B_M, IER)
224                    ELSE
225                       CALL CONV_DIF_PHI (X_S(1,M,LN), DIF_S(1,M,LN), &
226                         DISCRETIZE(7), U_S(1,M), V_S(1,M), W_S(1,M), &
227                         Flux_sSE, Flux_sSN, Flux_sST, M, A_M, B_M, IER)
228                    ENDIF
229     
230                    CALL BC_PHI (X_S(1,M,LN), BC_X_S(1,M,LN), &
231                       BC_XW_S(1,M,LN), BC_HW_X_S(1,M,LN), &
232                       BC_C_X_S(1,M,LN), M, A_M, B_M, IER)
233     
234                    CALL SOURCE_PHI (S_P, S_C, EPS, X_S(1,M,LN), M, A_M, B_M)
235     
236     ! Add point sources.
237                    IF(POINT_SOURCE) CALL POINT_SOURCE_PHI (X_S(1,M,LN), &
238                       PS_X_S(:,M,LN), PS_MASSFLOW_S(:,M), M, A_M, B_M, IER)
239     
240                    CALL CALC_RESID_S (X_S(1,M,LN), A_M, B_M, M, &
241                       NUM_RESID(RESID_X+(LN-1),M), &
242                       DEN_RESID(RESID_X+(LN-1),M), RESID(RESID_X+(LN-1),&
243                       M), MAX_RESID(RESID_X+(LN-1),M), IJK_RESID(RESID_X+(LN-1),M), &
244                       ZERO_X_GS, IER)
245     
246                    CALL UNDER_RELAX_S (X_S(1,M,LN), A_M, B_M, M, UR_FAC(7), IER)
247     
248     !               call check_ab_m(a_m, b_m, m, .false., ier)
249     !               write(*,*) resid(resid_x+(LN-1), m), &
250     !                  max_resid(resid_x+(LN-1), m), &
251     !                  ijk_resid(resid_x+(LN-1), m)
252     !               call write_ab_m(a_m, b_m, ijkmax2, m, ier)
253     !
254     !               call test_lin_eq(ijkmax2, ijmax2, imax2, a_m(1,-3,M),&
255     !                  1, DO_K, ier)
256     
257                    CALL ADJUST_LEQ (RESID(RESID_X+(LN-1),M), LEQ_IT(7), &
258                       LEQ_METHOD(7), LEQI, LEQM, IER)
259     
260                    CALL SOLVE_LIN_EQ ('X_s', 7, X_S(1,M,LN), A_M, B_M, M,&
261                       LEQI, LEQM, LEQ_SWEEP(7), LEQ_TOL(7), LEQ_PC(7), IER)
262     
263     ! Check for linear solver divergence.
264                    IF(ier == -2) Err_l(myPE) = 131
265     
266                    CALL BOUND_X (X_S(1,M,LN), IJKMAX2, IER)
267     !               call out_array(X_s(1,m,LN), 'X_s')
268     
269                 END DO
270     
271                 if(chi_scheme) call unset_chi()
272              ENDIF ! check for any species in phase m
273           END DO ! for m = 1, mmax
274     ! end solids phases species equations
275     ! ----------------------------------------------------------------<<<
276     
277           call unlock_ambm
278           call unlock_tmp_array
279     
280     
281     ! If the linear solver diverged, species mass fractions may take on
282     ! unphysical values. To prevent them from propogating through the domain
283     ! or causing failure in other routines, force an exit from iterate and
284     ! reduce the time step.
285           CALL global_all_sum(Err_l, Err_g)
286           IER = maxval(Err_g)
287     
288           RETURN
289           END SUBROUTINE SOLVE_SPECIES_EQ
290