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