File: N:\mfix\model\solve_scalar_eq.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: SOLVE_SCALAR_EQ                                         C
4     !  Purpose: Solve scalar transport equations                           C
5     !                                                                      C
6     !                                                                      C
7     !  Author: M. Syamlal                                 Date: 4-12-99    C
8     !                                                                      C
9     !                                                                      C
10     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
11           SUBROUTINE SOLVE_Scalar_EQ(IER)
12     
13     ! Modules
14     !---------------------------------------------------------------------//
15           use ambm, only: a_m, b_m, lock_ambm, unlock_ambm
16           use bc, only: bc_scalar, bc_hw_scalar, bc_c_scalar, bc_scalarw
17           use compar, only: ijkstart3, ijkend3
18           use fldvar, only: rop_g, rop_go, ep_g, u_g, v_g, w_g
19           use fldvar, only: rop_s, rop_so, ep_s, u_s, v_s, w_s
20           use fldvar, only: scalar, scalaro
21           use functions, only: fluid_at, zmax
22           use geometry, only: vol, ijkmax2
23           use leqsol, only: leq_it, leq_method, leq_sweep, leq_pc, leq_tol
24           use mflux, only: flux_ge, flux_gn, flux_gt
25           use mflux, only: flux_se, flux_sn, flux_st
26           use mflux, only: flux_gse, flux_gsn, flux_gst
27           use mflux, only: flux_sse, flux_ssn, flux_sst
28           use param, only: dimension_3
29           use param1, only: zero
30           use residual, only: resid, max_resid, ijk_resid
31           use residual, only: num_resid, den_resid
32           use residual, only: resid_sc
33           use run, only: discretize, added_mass, M_AM, odt
34           use rxns, only: sum_r_g, sum_r_s
35           use scalars, only: nscalar, phase4scalar
36           use scalars, only: scalar_c, scalar_p, dif_scalar
37           use toleranc, only: zero_ep_s
38           use ur_facs, only: ur_fac
39           use usr_src, only: call_usr_source, calc_usr_source
40           use usr_src, only: usr_scalar
41           IMPLICIT NONE
42     
43     ! Dummy arguments
44     !---------------------------------------------------------------------//
45     ! Error index
46           INTEGER, INTENT(INOUT) :: IER
47     
48     ! Local variables
49     !---------------------------------------------------------------------//
50     ! phase index
51           INTEGER :: M
52     ! species index
53           INTEGER :: NN
54     ! Indices
55           INTEGER :: IJK
56     !
57           DOUBLE PRECISION :: APO
58     !
59     ! linear equation solver method and iterations
60           INTEGER :: LEQM, LEQI
61     ! temporary variables in residual computation
62           DOUBLE PRECISION :: res1, mres1, num_res, den_res
63           INTEGER :: ires1
64     ! source vector: coefficient of dependent variable
65     ! becomes part of a_m matrix; must be positive
66           DOUBLE PRECISION :: S_P(DIMENSION_3)
67     ! source vector: constant part becomes part of b_m vector
68           DOUBLE PRECISION :: S_C(DIMENSION_3)
69     !
70           DOUBLE PRECISION :: EPS(DIMENSION_3)
71     !
72           character(LEN=8) :: Vname
73     !---------------------------------------------------------------------//
74     
75           call lock_ambm
76     
77           RESID(RESID_sc,0) = ZERO
78           NUM_RESID(RESID_sc,0) = ZERO
79           DEN_RESID(RESID_sc,0) = ZERO
80           MAX_RESID(RESID_sc,0) = ZERO
81           IJK_RESID(RESID_sc,0) = 0
82     
83     ! Loop over number of scalar equations
84           DO NN = 1, NScalar
85     
86     ! Setting M to the phase assigned to convecting the indicated scalar
87     ! equation
88              M = Phase4Scalar(NN)
89              CALL INIT_AB_M (A_M, B_M, IJKMAX2, M)
90     
91     ! Gas phase
92     ! ---------------------------------------------------------------->>>
93              IF(M == 0) THEN
94     
95                 DO IJK = IJKSTART3, IJKEND3
96                    IF (FLUID_AT(IJK)) THEN
97                       APO = ROP_GO(IJK)*VOL(IJK)*ODT
98                       S_P(IJK) = APO + (ZMAX(SUM_R_G(IJK)) + &
99                                              Scalar_p(IJK,NN))*VOL(IJK)
100                       S_C(IJK) = APO*ScalarO(IJK,NN) + &
101                                  Scalar(IJK,NN)*ZMAX((-SUM_R_G(IJK)))*VOL(IJK) + &
102                                  Scalar_c(IJK, NN)*VOL(IJK)
103                    ELSE
104                       S_P(IJK) = ZERO
105                       S_C(IJK) = ZERO
106                    ENDIF
107                 ENDDO
108     
109                 IF(.NOT.ADDED_MASS) THEN
110                    CALL CONV_DIF_PHI (Scalar(1,NN), DIF_Scalar(1,NN), &
111                                       DISCRETIZE(9), U_G, V_G, W_G, &
112                                       Flux_gE, Flux_gN, Flux_gT, M, A_M, B_M)
113                 ELSE
114                    CALL CONV_DIF_PHI (Scalar(1,NN), DIF_Scalar(1,NN), &
115                                       DISCRETIZE(9), U_G, V_G, W_G, &
116                                       Flux_gSE, Flux_gSN, Flux_gST, M, A_M, B_M)
117                 ENDIF
118     
119                 CALL BC_PHI (Scalar(1,NN), BC_Scalar(1,NN), &
120                              BC_ScalarW(1,NN), BC_HW_Scalar(1,NN), &
121                              BC_C_Scalar(1,NN), M, A_M, B_M)
122     
123                 CALL SOURCE_PHI (S_P, S_C, EP_G, Scalar(1,NN), M, A_M, B_M)
124     
125     ! calculate any usr source terms
126                 IF (CALL_USR_SOURCE(9)) CALL CALC_USR_SOURCE(USR_SCALAR,&
127                                       A_M, B_M, lM=M, lN=NN)
128     
129                 CALL CALC_RESID_S (Scalar(1,NN), A_M, B_M, M, num_res, &
130                                    den_res, res1, mres1, ires1, ZERO)
131                 RESID(RESID_sc,0) = RESID(RESID_sc,0)+res1
132                 NUM_RESID(RESID_sc,0) = NUM_RESID(RESID_sc,0)+num_res
133                 DEN_RESID(RESID_sc,0) = DEN_RESID(RESID_sc,0)+den_res
134                 if(mres1 .gt. MAX_RESID(RESID_sc,0))then
135                    MAX_RESID(RESID_sc,0) = mres1
136                    IJK_RESID(RESID_sc,0) = ires1
137                 endif
138     
139                 CALL UNDER_RELAX_S (Scalar(1,NN), A_M, B_M, M, UR_FAC(9))
140     
141     !           call check_ab_m(a_m, b_m, m, .false., ier)
142     !           call write_ab_m(a_m, b_m, ijkmax2, m, ier)
143     !           write(*,*) res1, mres1, ires1
144     !           call test_lin_eq(ijkmax2, ijmax2, imax2, a_m(1, -3, 0), 1, &
145     !                            DO_K, ier)
146     
147                 CALL ADJUST_LEQ (res1, LEQ_IT(9), LEQ_METHOD(9), &
148                    LEQI, LEQM)
149     
150                 write(Vname, '(A,I2)')'Scalar',NN
151                 CALL SOLVE_LIN_EQ (Vname, 9, Scalar(1,NN), A_M, B_M, M,&
152                                    LEQI, LEQM, LEQ_SWEEP(9), LEQ_TOL(9),&
153                                    LEQ_PC(9), IER)
154     !           call out_array(Scalar(1, N), Vname)
155     
156     ! Solids phase
157     ! ---------------------------------------------------------------->>>
158              ELSE   ! (m/=0) , i.e., solids phase
159     
160                 DO IJK = IJKSTART3, IJKEND3
161                    IF (FLUID_AT(IJK)) THEN
162                       APO = ROP_sO(IJK, M)*VOL(IJK)*ODT
163                       S_P(IJK) = APO + (ZMAX(SUM_R_s(IJK, M)) + &
164                                              Scalar_p(IJK, NN))*VOL(IJK)
165                       S_C(IJK) = APO*ScalarO(IJK,NN) + &
166                                  Scalar(IJK,NN)*ZMAX((-SUM_R_s(IJK, M)))*VOL(IJK)+&
167                                  Scalar_c(IJK, NN)*VOL(IJK)
168                       EPs(IJK) = EP_s(IJK, M)
169                    ELSE
170                       S_P(IJK) = ZERO
171                       S_C(IJK) = ZERO
172                       EPS(IJK) = ZERO
173                    ENDIF
174                 ENDDO
175     
176                 IF(.NOT.ADDED_MASS .OR. M /= M_AM) THEN
177                    CALL CONV_DIF_PHI (Scalar(1,NN), DIF_Scalar(1,NN), &
178                                       DISCRETIZE(9), U_s(1,m), V_s(1,m), &
179                                       W_s(1,m), Flux_sE(1,M), Flux_sN(1,M),&
180                                       Flux_sT(1,M), M, A_M, B_M)
181                 ELSE ! virtual mass term added for M = M_AM ONLY!!!!
182                    CALL CONV_DIF_PHI (Scalar(1,NN), DIF_Scalar(1,NN), &
183                                       DISCRETIZE(9), U_s(1,m), V_s(1,m), &
184                                       W_s(1,m), Flux_sSE, Flux_sSN, Flux_sST, &
185                                       M, A_M, B_M)
186                 ENDIF
187     
188                 CALL BC_PHI (Scalar(1,NN), BC_Scalar(1,NN), BC_ScalarW(1,NN), &
189                              BC_HW_Scalar(1,NN), BC_C_Scalar(1,NN), M, A_M, B_M)
190     
191                 CALL SOURCE_PHI (S_P, S_C, EPs, Scalar(1,NN), M, A_M, B_M)
192     
193     ! calculate any usr source terms
194                 IF (CALL_USR_SOURCE(9)) CALL CALC_USR_SOURCE(USR_SCALAR,&
195                                       A_M, B_M, lM=M, lN=NN)
196     
197     
198                 CALL CALC_RESID_S (Scalar(1,NN), A_M, B_M, M, num_res, &
199                                   den_res, res1, mres1, ires1, ZERO)
200                 RESID(RESID_sc,0) = RESID(RESID_sc,0)+res1
201                 NUM_RESID(RESID_sc,0) = NUM_RESID(RESID_sc,0)+num_res
202                 DEN_RESID(RESID_sc,0) = DEN_RESID(RESID_sc,0)+den_res
203                 if(mres1 .gt. MAX_RESID(RESID_sc,0))then
204                    MAX_RESID(RESID_sc,0) = mres1
205                    IJK_RESID(RESID_sc,0) = ires1
206                 endif
207     
208                 CALL UNDER_RELAX_S (Scalar(1,NN), A_M, B_M, M, UR_FAC(9))
209     
210     !            call check_ab_m(a_m, b_m, m, .false., ier)
211     !            call write_ab_m(a_m, b_m, ijkmax2, m, ier)
212     !            write(*,*) res1, mres1, ires1
213     !            call test_lin_eq(ijkmax2, ijmax2, imax2, a_m(1, -3, 0), 1, &
214     !                             DO_K, ier)
215     
216                 CALL ADJUST_LEQ (res1, LEQ_IT(9), LEQ_METHOD(9), &
217                                  LEQI, LEQM)
218     
219                 write(Vname, '(A,I2)')'Scalar',NN
220                 CALL SOLVE_LIN_EQ (Vname, 9, Scalar(1,NN), A_M, B_M, M, &
221                                    LEQI, LEQM, LEQ_SWEEP(9), LEQ_TOL(9), &
222                                    LEQ_PC(9), IER)
223     !            call out_array(Scalar(1, N), Vname)
224     
225              ENDIF
226           ENDDO
227     
228           call unlock_ambm
229     
230           RETURN
231           END SUBROUTINE SOLVE_Scalar_EQ
232