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