1 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC 2 ! C 3 ! Subroutine: CONV_ROP_s C 4 ! Purpose: Determine convection terms for solids continuity C 5 ! equation. Master routine. C 6 ! C 7 ! Author: M. Syamlal Date: 18-MAR-97 C 8 ! Reviewer: Date: C 9 ! C 10 ! Literature/Document References: C 11 ! Variables referenced: C 12 ! Variables modified: C 13 ! Local variables: C 14 ! C 15 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C 16 17 SUBROUTINE CONV_ROP_S(A_M, B_M, M, IER) 18 19 !----------------------------------------------- 20 ! Modules 21 !----------------------------------------------- 22 USE param 23 USE param1 24 USE fldvar 25 USE run 26 USE compar 27 IMPLICIT NONE 28 !----------------------------------------------- 29 ! Dummy arguments 30 !----------------------------------------------- 31 ! Septadiagonal matrix A_m 32 DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M) 33 ! Vector b_m 34 DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M) 35 ! Solids phase index 36 INTEGER, INTENT(IN) :: M 37 ! Error index 38 INTEGER, INTENT(INOUT) :: IER 39 !----------------------------------------------- 40 41 IF (DISCRETIZE(2) == 0) THEN ! 0 & 1 => first order upwinding 42 CALL CONV_ROP_S0 (A_M, M) 43 ELSE 44 CALL CONV_ROP_S1 (A_M, M) 45 ENDIF 46 47 RETURN 48 END SUBROUTINE CONV_ROP_S 49 50 51 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC 52 ! C 53 ! Subroutine: CONV_ROP_s0 C 54 ! Purpose: Determine convection terms for solids continuity equation. C 55 ! C 56 ! Notes: The off-diagonal coefficients calculated here must be C 57 ! positive. The center coefficient and the source vector C 58 ! are negative. -- First order upwinding C 59 ! C 60 ! Author: M. Syamlal Date: 3-JUL-96 C 61 ! Reviewer: Date: C 62 ! C 63 ! Literature/Document References: C 64 ! Variables referenced: C 65 ! Variables modified: C 66 ! Local variables: C 67 ! C 68 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C 69 70 SUBROUTINE CONV_ROP_S0(A_M, M) 71 72 !----------------------------------------------- 73 ! Modules 74 !----------------------------------------------- 75 USE param 76 USE param1 77 USE fldvar 78 USE run 79 USE parallel 80 USE matrix 81 USE physprop 82 USE geometry 83 USE indices 84 USE pgcor 85 USE pscor 86 USE compar 87 USE functions 88 IMPLICIT NONE 89 !----------------------------------------------- 90 ! Dummy arguments 91 !----------------------------------------------- 92 ! Septadiagonal matrix A_m 93 DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M) 94 ! Solids phase index 95 INTEGER, INTENT(IN) :: M 96 !----------------------------------------------- 97 ! Local variables 98 !----------------------------------------------- 99 ! Indices 100 INTEGER :: I, J, K, IJK, IPJK, IJPK, IJKP 101 INTEGER :: IMJK, IJMK, IJKM 102 !----------------------------------------------- 103 104 ! Calculate convection-diffusion fluxes through each of the faces 105 106 !!$omp parallel do private( I, J, K, IJK, IPJK, IJPK, IJKP, & 107 !!$omp& IMJK, IJMK, IJKM) & 108 !!$omp& schedule(static) 109 DO IJK = ijkstart3, ijkend3 110 IF (PHASE_4_P_G(IJK)/=M .AND. PHASE_4_P_S(IJK)/=M) THEN 111 ! if either phase_4_p_g or phase_4_p_s are assigned the index of the 112 ! current solids phase then this section is skipped. currently, 113 ! phase_4_p_g and phase_4_p_s are always assigned to 0 and undefined, 114 ! respectively. Hence this section should always be entered. 115 116 ! it was previously possible for phase_4_p_g to be assigned the index of 117 ! the gas phase in some cells, and the index of a solids phase in other 118 ! cells if that solids phase has close_packed=F and was in higher 119 ! concentrations than the gas. in such a case this branch would become 120 ! skipped, while the corresponding gas continuity would become 121 ! activated. moreover, that solids phase's volume fraction would be 122 ! corrected rather than the gas phases void fraction in calc_vol_fr. 123 124 ! if phase_4_p_s were to be assigned the index of a solids phase that 125 ! can close pack if the current cell exhibited close packing, then 126 ! this branch would become skipped. moreover, that solids phase's 127 ! volume fraction would then be corected based on the value of 128 ! maximum packing in that cell and the sum of all other solids that 129 ! can close pack in calc_vol_fr. 130 I = I_OF(IJK) 131 J = J_OF(IJK) 132 K = K_OF(IJK) 133 IPJK = IP_OF(IJK) 134 IJPK = JP_OF(IJK) 135 IJKP = KP_OF(IJK) 136 IMJK = IM_OF(IJK) 137 IJMK = JM_OF(IJK) 138 IJKM = KM_OF(IJK) 139 140 ! East face (i+1/2, j, k) 141 A_M(IJK,E,M) = ZMAX((-U_S(IJK,M)))*AYZ(IJK) 142 A_M(IPJK,W,M) = ZMAX(U_S(IJK,M))*AYZ(IJK) 143 144 ! North face (i, j+1/2, k) 145 A_M(IJK,N,M) = ZMAX((-V_S(IJK,M)))*AXZ(IJK) 146 A_M(IJPK,S,M) = ZMAX(V_S(IJK,M))*AXZ(IJK) 147 148 ! Top face (i, j, k+1/2) 149 IF (DO_K) THEN 150 A_M(IJK,T,M) = ZMAX((-W_S(IJK,M)))*AXY(IJK) 151 A_M(IJKP,B,M) = ZMAX(W_S(IJK,M))*AXY(IJK) 152 ENDIF 153 154 ! Modify west (i-1/2,j,k), south (i j-1/2,k) and bottom (i,j,k-1/2) 155 ! faces if the neighboring west, south, bottom cells have 156 ! phase_4_p_g of m or phase_4_p_s of m. 157 IF (PHASE_4_P_G(IMJK)==M .OR. PHASE_4_P_S(IMJK)==M) & 158 A_M(IJK,W,M) = ZMAX(U_S(IMJK,M))*AYZ(IMJK) 159 IF (PHASE_4_P_G(IJMK)==M .OR. PHASE_4_P_S(IJMK)==M) & 160 A_M(IJK,S,M) = ZMAX(V_S(IJMK,M))*AXZ(IJMK) 161 162 IF (DO_K) THEN 163 IF (PHASE_4_P_G(IJKM)==M .OR. PHASE_4_P_S(IJKM)==M) & 164 A_M(IJK,B,M) = ZMAX(W_S(IJKM,M))*AXY(IJKM) 165 ENDIF 166 ENDIF 167 END DO 168 169 RETURN 170 END SUBROUTINE CONV_ROP_S0 171 172 173 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC 174 ! C 175 ! Subroutine: CONV_ROP_s1 C 176 ! Purpose: Determine convection terms for solids continuity equation. C 177 ! Notes: The off-diagonal coefficients calculated here must be C 178 ! positive. The center coefficient and the source vector C 179 ! are negative. -- Higher order methods C 180 ! See conv_rop_s0 for additional details. C 181 ! C 182 ! Author: M. Syamlal Date: 18-MAR-97 C 183 ! Reviewer: Date: C 184 ! C 185 ! C 186 ! Literature/Document References: C 187 ! C 188 ! Variables referenced: C 189 ! Variables modified: C 190 ! C 191 ! Local variables: C 192 ! C 193 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C 194 195 SUBROUTINE CONV_ROP_S1(A_M, M) 196 197 !----------------------------------------------- 198 ! Modules 199 !----------------------------------------------- 200 USE param 201 USE param1 202 USE fldvar 203 USE run 204 USE parallel 205 USE matrix 206 USE physprop 207 USE geometry 208 USE indices 209 USE pgcor 210 USE pscor 211 Use xsi_array 212 USE compar 213 USE xsi 214 USE functions 215 IMPLICIT NONE 216 !----------------------------------------------- 217 ! Dummy arguments 218 !----------------------------------------------- 219 ! Septadiagonal matrix A_m 220 DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M) 221 ! Solids phase index 222 INTEGER, INTENT(IN) :: M 223 !----------------------------------------------- 224 ! Local variables 225 !----------------------------------------------- 226 ! Indices 227 INTEGER :: I, J, K, IJK, IPJK, IJPK, IJKP 228 INTEGER :: IMJK, IJMK, IJKM 229 ! loezos 230 INTEGER :: incr 231 ! temporary use of global arrays: 232 ! xsi_array: convection weighting factors 233 ! DOUBLE PRECISION :: XSI_e(DIMENSION_3), XSI_n(DIMENSION_3),& 234 ! XSI_t(DIMENSION_3) 235 !----------------------------------------------- 236 237 call lock_xsi_array 238 239 ! Loezos: 240 incr=0 241 242 ! Calculate convection factors 243 CALL CALC_XSI(DISCRETIZE(2), ROP_S(1,M), U_S(1,M), V_S(1,M),& 244 W_S(1,M), XSI_E, XSI_N, XSI_T, incr) 245 246 ! Calculate convection-diffusion fluxes through each of the faces 247 248 !!$omp parallel do private( I, J, K, IJK, IPJK, IJPK, IJKP, & 249 !!$omp& IMJK, IJMK, IJKM) & 250 !!$omp& schedule(static) 251 DO IJK = ijkstart3, ijkend3 252 IF (PHASE_4_P_G(IJK)/=M .AND. PHASE_4_P_S(IJK)/=M) THEN 253 I = I_OF(IJK) 254 J = J_OF(IJK) 255 K = K_OF(IJK) 256 IPJK = IP_OF(IJK) 257 IJPK = JP_OF(IJK) 258 IJKP = KP_OF(IJK) 259 IMJK = IM_OF(IJK) 260 IJMK = JM_OF(IJK) 261 IJKM = KM_OF(IJK) 262 263 ! East face (i+1/2, j, k) 264 A_M(IJK,E,M) = -XSI_E(IJK)*U_S(IJK,M)*AYZ(IJK) 265 A_M(IPJK,W,M) = (ONE - XSI_E(IJK))*U_S(IJK,M)*AYZ(IJK) 266 267 ! North face (i, j+1/2, k) 268 A_M(IJK,N,M) = -XSI_N(IJK)*V_S(IJK,M)*AXZ(IJK) 269 A_M(IJPK,S,M) = (ONE - XSI_N(IJK))*V_S(IJK,M)*AXZ(IJK) 270 271 ! Top face (i, j, k+1/2) 272 IF (DO_K) THEN 273 A_M(IJK,T,M) = -XSI_T(IJK)*W_S(IJK,M)*AXY(IJK) 274 A_M(IJKP,B,M) = (ONE - XSI_T(IJK))*W_S(IJK,M)*AXY(IJK) 275 ENDIF 276 277 IF (PHASE_4_P_G(IMJK)==M .OR. PHASE_4_P_S(IMJK)==M) A_M(IJK,W,M) = & 278 (ONE - XSI_E(IMJK))*U_S(IMJK,M)*AYZ(IMJK) 279 IF (PHASE_4_P_G(IJMK)==M .OR. PHASE_4_P_S(IJMK)==M) A_M(IJK,S,M) = & 280 (ONE - XSI_N(IJMK))*V_S(IJMK,M)*AXZ(IJMK) 281 IF (DO_K) THEN 282 IF (PHASE_4_P_G(IJKM)==M .OR. PHASE_4_P_S(IJKM)==M) A_M(IJK,B,M)& 283 = (ONE - XSI_T(IJKM))*W_S(IJKM,M)*AXY(IJKM) 284 ENDIF 285 ENDIF 286 287 ENDDO 288 289 call unlock_xsi_array 290 291 RETURN 292 END SUBROUTINE CONV_ROP_S1 293 294