File: RELATIVE:/../../../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 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