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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Module name: SOLVE_SCALAR_EQ(IER)                                   C
4     !  Purpose: Solve scalar transport equations                           C
5     !                                                                      C
6     !                                                                      C
7     !  Author: M. Syamlal                                 Date: 4-12-99    C
8     !  Reviewer:                                          Date:            C
9     !                                                                      C
10     !                                                                      C
11     !  Literature/Document References:                                     C
12     !                                                                      C
13     !  Variables referenced:                                               C
14     !  Variables modified:                                                 C
15     !                                                                      C
16     !  Local variables:                                                    C
17     !                                                                      C
18     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
19     !
20           SUBROUTINE SOLVE_Scalar_EQ(IER)
21     !...Translated by Pacific-Sierra Research VAST-90 2.06G5  12:17:31  12/09/98
22     !...Switches: -xf
23     !
24     !  Include param.inc file to specify parameter values
25     !
26     !-----------------------------------------------
27     !   M o d u l e s
28     !-----------------------------------------------
29           USE param
30           USE param1
31           USE toleranc
32           USE run
33           USE physprop
34           USE geometry
35           USE fldvar
36           USE output
37           USE indices
38           USE drag
39           USE residual
40           USE ur_facs
41           USE pgcor
42           USE pscor
43           USE leqsol
44           USE bc
45           USE energy
46           USE rxns
47           USE scalars
48           Use ambm
49           Use tmp_array, S_p => Array1, S_c => Array2, EPs => Array3, VxGama => Array4
50           USE compar
51           USE mflux
52           USE functions
53           IMPLICIT NONE
54     !-----------------------------------------------
55     !   G l o b a l   P a r a m e t e r s
56     !-----------------------------------------------
57     !-----------------------------------------------
58     !   D u m m y   A r g u m e n t s
59     !-----------------------------------------------
60     !
61     !                      Error index
62           INTEGER          IER
63     !-----------------------------------------------
64     !   L o c a l   P a r a m e t e r s
65     !-----------------------------------------------
66     !-----------------------------------------------
67     !   L o c a l   V a r i a b l e s
68     !-----------------------------------------------
69     !
70     !
71     !                      phase index
72           INTEGER          m
73     !
74     !                      species index
75           INTEGER          n
76     !
77           DOUBLE PRECISION apo
78     !
79     
80     !
81     !                      temporary variables in residual computation
82           DOUBLE PRECISION res1, mres1, num_res, den_res
83           INTEGER          ires1
84     !
85     !                      Indices
86           INTEGER          IJK
87     !
88     !                      linear equation solver method and iterations
89           INTEGER          LEQM, LEQI
90     
91           character(LEN=8) :: Vname
92     
93     !-----------------------------------------------
94     
95           call lock_ambm
96           call lock_tmp_array
97     
98           RESID(RESID_sc,0) = ZERO
99           NUM_RESID(RESID_sc,0) = ZERO
100           DEN_RESID(RESID_sc,0) = ZERO
101           MAX_RESID(RESID_sc,0) = ZERO
102           IJK_RESID(RESID_sc,0) = 0
103     !
104     !     Fluid phase species mass balance equations
105     !
106            DO N = 1, NScalar
107     
108               M = Phase4Scalar(N)
109     
110               CALL INIT_AB_M (A_M, B_M, IJKMAX2, M, IER)
111     
112               IF(M == 0) THEN
113     
114               DO IJK = IJKSTART3, IJKEND3
115     !
116                    IF (FLUID_AT(IJK)) THEN
117                       APO = ROP_GO(IJK)*VOL(IJK)*ODT
118                       S_P(IJK) = APO + (  ZMAX(SUM_R_G(IJK)) &
119                                         + Scalar_p(IJK, N)      )*VOL(IJK)
120                       S_C(IJK) =   APO*ScalarO(IJK,N) &
121                                  + Scalar(IJK,N)*ZMAX((-SUM_R_G(IJK)))*VOL(IJK) &
122                                  + Scalar_c(IJK, N)*VOL(IJK)
123                    ELSE
124     !
125                       S_P(IJK) = ZERO
126                       S_C(IJK) = ZERO
127     !
128                    ENDIF
129                 END DO
130     
131     
132                 IF(.NOT.ADDED_MASS) THEN
133                    CALL CONV_DIF_PHI (Scalar(1,N), DIF_Scalar(1,N), DISCRETIZE(9), &
134                                    U_G, V_G, W_G, Flux_gE, Flux_gN, Flux_gT, M, A_M, B_M, IER)
135                 ELSE
136                    CALL CONV_DIF_PHI (Scalar(1,N), DIF_Scalar(1,N), DISCRETIZE(9), &
137                                    U_G, V_G, W_G, Flux_gSE, Flux_gSN, Flux_gST, M, A_M, B_M, IER)
138                 ENDIF
139     !
140     !
141                 CALL BC_PHI (Scalar(1,N), BC_Scalar(1,N), BC_ScalarW(1,N), BC_HW_Scalar(1,N), &
142                              BC_C_Scalar(1,N), M, A_M, B_M, IER)
143     !
144     !
145                 CALL SOURCE_PHI (S_P, S_C, EP_G, Scalar(1,N), M, A_M, B_M, IER)
146     !
147                 CALL CALC_RESID_S (Scalar(1,N), A_M, B_M, M, num_res, den_res, res1, &
148                    mres1, ires1, ZERO, IER)
149                    RESID(RESID_sc,0) = RESID(RESID_sc,0)+res1
150                    NUM_RESID(RESID_sc,0) = NUM_RESID(RESID_sc,0)+num_res
151                    DEN_RESID(RESID_sc,0) = DEN_RESID(RESID_sc,0)+den_res
152                    if(mres1 .gt. MAX_RESID(RESID_sc,0))then
153                      MAX_RESID(RESID_sc,0) = mres1
154                      IJK_RESID(RESID_sc,0) = ires1
155                    endif
156     !
157                 CALL UNDER_RELAX_S (Scalar(1,N), A_M, B_M, M, UR_FAC(9), IER)
158     !
159     !          call check_ab_m(a_m, b_m, m, .false., ier)
160     !          call write_ab_m(a_m, b_m, ijkmax2, m, ier)
161     !          write(*,*) res1, mres1, &
162     !           ires1
163     !
164     !          call test_lin_eq(ijkmax2, ijmax2, imax2, a_m(1, -3, 0), 1, DO_K, &
165     !          ier)
166     !
167                 CALL ADJUST_LEQ (res1, LEQ_IT(9), LEQ_METHOD(9), &
168                    LEQI, LEQM, IER)
169     !
170                 write(Vname, '(A,I2)')'Scalar',N
171                 CALL SOLVE_LIN_EQ (Vname, 9, Scalar(1,N), A_M, B_M, M, LEQI, LEQM, &
172                                  LEQ_SWEEP(9), LEQ_TOL(9), LEQ_PC(9), IER)
173     !          call out_array(Scalar(1, N), Vname)
174     !
175              ELSE
176     
177                 DO IJK = IJKSTART3, IJKEND3
178     !
179                    IF (FLUID_AT(IJK)) THEN
180                       APO = ROP_sO(IJK, M)*VOL(IJK)*ODT
181                       S_P(IJK) = APO + (  ZMAX(SUM_R_s(IJK, M)) &
182                                         + Scalar_p(IJK, N)       )*VOL(IJK)
183                       S_C(IJK) =   APO*ScalarO(IJK,N) &
184                                  + Scalar(IJK,N)*ZMAX((-SUM_R_s(IJK, M)))*VOL(IJK)&
185                                  + Scalar_c(IJK, N)*VOL(IJK)
186                       EPs(IJK) = EP_s(IJK, M)
187                    ELSE
188     !
189                       S_P(IJK) = ZERO
190                       S_C(IJK) = ZERO
191                       EPS(IJK) = ZERO
192     !
193                    ENDIF
194                 END DO
195     
196     
197                 IF(.NOT.ADDED_MASS .OR. M /= M_AM) THEN
198                    CALL CONV_DIF_PHI (Scalar(1,N), DIF_Scalar(1,N), DISCRETIZE(9), &
199                                    U_s(1,m), V_s(1,m), W_s(1,m), Flux_sE(1,M), Flux_sN(1,M), Flux_sT(1,M), M, &
200                                    A_M, B_M, IER)
201                 ELSE ! virtual mass term added for M = M_AM ONLY!!!!
202                    CALL CONV_DIF_PHI (Scalar(1,N), DIF_Scalar(1,N), DISCRETIZE(9), &
203                                    U_s(1,m), V_s(1,m), W_s(1,m), Flux_sSE, Flux_sSN, Flux_sST, M, &
204                                    A_M, B_M, IER)
205                 ENDIF
206     !
207     !
208                 CALL BC_PHI (Scalar(1,N), BC_Scalar(1,N), BC_ScalarW(1,N), BC_HW_Scalar(1,N), &
209                              BC_C_Scalar(1,N), M, A_M, B_M, IER)
210     !
211     !
212                 CALL SOURCE_PHI (S_P, S_C, EPs, Scalar(1,N), M, A_M, B_M, IER)
213     !
214                 CALL CALC_RESID_S (Scalar(1,N), A_M, B_M, M, num_res, den_res, res1, &
215                    mres1, ires1, ZERO, IER)
216                    RESID(RESID_sc,0) = RESID(RESID_sc,0)+res1
217                    NUM_RESID(RESID_sc,0) = NUM_RESID(RESID_sc,0)+num_res
218                    DEN_RESID(RESID_sc,0) = DEN_RESID(RESID_sc,0)+den_res
219                    if(mres1 .gt. MAX_RESID(RESID_sc,0))then
220                      MAX_RESID(RESID_sc,0) = mres1
221                      IJK_RESID(RESID_sc,0) = ires1
222                    endif
223     !
224                 CALL UNDER_RELAX_S (Scalar(1,N), A_M, B_M, M, UR_FAC(9), IER)
225     !
226     !          call check_ab_m(a_m, b_m, m, .false., ier)
227     !          call write_ab_m(a_m, b_m, ijkmax2, m, ier)
228     !          write(*,*) res1, mres1, ires1
229     !
230     !          call test_lin_eq(ijkmax2, ijmax2, imax2, a_m(1, -3, 0), 1, DO_K, &
231     !          ier)
232     !
233                 CALL ADJUST_LEQ (res1, LEQ_IT(9), LEQ_METHOD(9), &
234                    LEQI, LEQM, IER)
235     !
236                 write(Vname, '(A,I2)')'Scalar',N
237                 CALL SOLVE_LIN_EQ (Vname, 9, Scalar(1,N), A_M, B_M, M, LEQI, LEQM, &
238                                  LEQ_SWEEP(9), LEQ_TOL(9), LEQ_PC(9), IER)
239     !          call out_array(Scalar(1, N), Vname)
240     !
241              END IF
242           END DO
243     
244           call unlock_ambm
245           call unlock_tmp_array
246     
247           RETURN
248           END SUBROUTINE SOLVE_Scalar_EQ
249     
250     
251     !// Comments on the modifications for DMP version implementation
252     !// 001 Include header file and common declarations for parallelization
253     !// 350 Changed do loop limits: 1,ijkmax2-> ijkstart3, ijkend3
254