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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: SOURCE_ROP_s                                            C
4     !  Purpose: Determine source terms for solids continuity equation.     C
5     !                                                                      C
6     !  Notes: The off-diagonal coefficients are positive. The center       C
7     !         coefficient and the source vector are negative.              C
8     !         See conv_rop_s0 for additional details.                      C
9     !                                                                      C
10     !  Author: M. Syamlal                                 Date: 3 -JUL-96  C
11     !  Reviewer:                                          Date:            C
12     !                                                                      C
13     !  Literature/Document References:                                     C
14     !  Variables referenced:                                               C
15     !  Variables modified:                                                 C
16     !  Local variables:                                                    C
17     !                                                                      C
18     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
19     
20           SUBROUTINE SOURCE_ROP_S(A_M, B_M, M)
21     
22     !-----------------------------------------------
23     ! Modules
24     !-----------------------------------------------
25           USE param
26           USE param1
27           USE parallel
28           USE matrix
29           USE fldvar
30           USE rxns
31           USE run
32           USE geometry
33           USE indices
34           USE pgcor
35           USE pscor
36           USE compar
37           USE functions
38           IMPLICIT NONE
39     !-----------------------------------------------
40     ! Dummy arguments
41     !-----------------------------------------------
42     ! Septadiagonal matrix A_m
43           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
44     ! Vector b_m
45           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
46     ! solids phase index
47           INTEGER, INTENT(IN) :: M
48     !-----------------------------------------------
49     ! Local variables
50     !-----------------------------------------------
51     ! DEL dot V
52           DOUBLE PRECISION :: DEL_V
53     ! Mass source
54           DOUBLE PRECISION :: Src
55     ! Indices
56           INTEGER :: I, J, K, IJK, IMJK, IJMK, IJKM
57     ! error message
58           CHARACTER(LEN=80) :: LINE
59     !-----------------------------------------------
60     
61     !!$omp  parallel do private( I, J, K, IJK, IMJK, IJMK, IJKM,  DEL_V, &
62     !!$omp&  Src, LINE) &
63     !!$omp&  schedule(static)
64           DO IJK = ijkstart3, ijkend3
65              IF (FLUID_AT(IJK) .AND. PHASE_4_P_G(IJK)/=M .AND. &
66                  PHASE_4_P_S(IJK)/=M) THEN
67     
68                 I = I_OF(IJK)
69                 J = J_OF(IJK)
70                 K = K_OF(IJK)
71                 IMJK = IM_OF(IJK)
72                 IJMK = JM_OF(IJK)
73                 IJKM = KM_OF(IJK)
74     
75                 DEL_V = U_S(IJK,M)*AYZ(IJK) - U_S(IMJK,M)*AYZ(IMJK) +&
76                         V_S(IJK,M)*AXZ(IJK) - V_S(IJMK,M)*AXZ(IJMK) + &
77                         W_S(IJK,M)*AXY(IJK) - W_S(IJKM,M)*AXY(IJKM)
78     
79                 IF (ROP_S(IJK,M) > ZERO) THEN
80                    SRC = VOL(IJK)*ZMAX((-SUM_R_S(IJK,M)))/ROP_S(IJK,M)
81                 ELSE
82                    SRC = ZERO
83                 ENDIF
84     
85                 A_M(IJK,0,M) = -(A_M(IJK,E,M)+A_M(IJK,W,M)+A_M(IJK,N,M)+&
86                                  A_M(IJK,S,M)+A_M(IJK,T,M)+A_M(IJK,B,M)+&
87                                  VOL(IJK)*ODT+ZMAX(DEL_V)+SRC)
88                 B_M(IJK,M) = -(ROP_SO(IJK,M)*VOL(IJK)*ODT+&
89                                ZMAX((-DEL_V))*ROP_S(IJK,M)+&
90                                ZMAX(SUM_R_S(IJK,M))*VOL(IJK))
91     
92                 IF (ABS(A_M(IJK,0,M)) < SMALL_NUMBER) THEN
93                    IF (ABS(B_M(IJK,M)) < SMALL_NUMBER) THEN
94                       A_M(IJK,0,M) = -ONE            ! Equation is undefined.
95                       B_M(IJK,M) = -ROP_S(IJK,M)     ! Use existing value
96                    ELSE
97     !!$omp             critical
98                       WRITE (LINE, '(A,I6,A,I1,A,G12.5)') 'Error: At IJK = ', IJK, &
99                          ' M = ', M, ' A = 0 and b = ', B_M(IJK,M)
100                       CALL WRITE_ERROR ('SOURCE_ROP_s', LINE, 1)
101     !!$omp             end critical
102                    ENDIF
103                 ENDIF
104              ELSE
105     ! set the value of rop_s in all wall and flow boundary cells to what is
106     ! known for that cell
107                 A_M(IJK,E,M) = ZERO
108                 A_M(IJK,W,M) = ZERO
109                 A_M(IJK,N,M) = ZERO
110                 A_M(IJK,S,M) = ZERO
111                 A_M(IJK,T,M) = ZERO
112                 A_M(IJK,B,M) = ZERO
113                 A_M(IJK,0,M) = -ONE
114                 B_M(IJK,M) = -ROP_S(IJK,M)
115              ENDIF
116           ENDDO
117     
118           RETURN
119           END SUBROUTINE SOURCE_ROP_S
120     
121     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
122     !                                                                      C
123     !  Subroutine: POINT_SOURCE_ROP_S                                      C
124     !  Purpose: Adds point sources to the solids continuity equation.      C
125     !                                                                      C
126     !                                                                      C
127     !  Author: J. Musser                                  Date: 10-JUN-13  C
128     !  Reviewer:                                          Date:            C
129     !                                                                      C
130     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
131           SUBROUTINE POINT_SOURCE_ROP_S(B_M, M)
132     
133           use compar
134           use constant
135           use geometry
136           use indices
137           use physprop
138           use ps
139           use run
140           use functions
141     
142           IMPLICIT NONE
143     !-----------------------------------------------
144     ! Dummy arguments
145     !-----------------------------------------------
146     ! Vector b_m
147           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
148     ! Solids phase index.
149           INTEGER, INTENT(IN) :: M
150     
151     !-----------------------------------------------
152     ! Local Variables
153     !-----------------------------------------------
154     
155     ! Indices
156           INTEGER :: IJK, I, J, K
157           INTEGER :: PSV
158     
159     ! terms of bm expression
160           DOUBLE PRECISION :: pSource
161     
162     !-----------------------------------------------
163           PS_LP: do PSV = 1, DIMENSION_PS
164     
165              if(.NOT.PS_DEFINED(PSV)) cycle PS_LP
166              if(PS_MASSFLOW_S(PSV,M) < small_number) cycle PS_LP
167     
168              do k = PS_K_B(PSV), PS_K_T(PSV)
169              do j = PS_J_S(PSV), PS_J_N(PSV)
170              do i = PS_I_W(PSV), PS_I_E(PSV)
171     
172                 if(.NOT.IS_ON_myPE_plus2layers(I,J,K)) cycle
173                 IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
174     
175                 ijk = funijk(i,j,k)
176                 if(fluid_at(ijk)) then
177     
178                    pSource =  PS_MASSFLOW_S(PSV,M) * &
179                       (VOL(IJK)/PS_VOLUME(PSV))
180     
181                    B_M(IJK,M) = B_M(IJK,M) - pSource
182                 endif
183              enddo
184              enddo
185              enddo
186     
187           enddo PS_LP
188     
189           RETURN
190           END SUBROUTINE POINT_SOURCE_ROP_S
191