File: RELATIVE:/../../../mfix.git/model/des/dif_phi_source_des.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !                                                                      !
3     !  Subroutine: SOURCE_phi                                              !
4     !  Author: M. Syamlal                                 Date: 30-APR-97  !
5     !                                                                      !
6     !  Purpose: Determine source terms for phi eq. The terms               !
7     !     appear in the center coefficient and RHS vector. The center      !
8     !     coefficient and source vector are negative. The off-diagonal     !
9     !     coefficients are positive. S_p must be positive.                 !
10     !                                                                      !
11     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
12           SUBROUTINE DIF_PHI_SOURCE_DES(PHI, M, A_M, B_M, lDT)
13     
14     !-----------------------------------------------
15     ! Modules
16     !-----------------------------------------------
17           USE param
18           USE param1
19           USE parallel
20           USE matrix
21           USE scales
22           USE physprop
23           USE fldvar
24           USE rxns
25           USE run
26           USE toleranc
27           USE geometry
28           USE indices
29           USE compar
30           USE fun_avg
31           USE functions
32     
33           IMPLICIT NONE
34     !-----------------------------------------------
35     ! Dummy Arguments
36     !-----------------------------------------------
37     ! Source term on RHS
38           DOUBLE PRECISION :: S_C
39     ! Phi
40           DOUBLE PRECISION, INTENT(IN) :: Phi(DIMENSION_3)
41     ! phase index
42           INTEGER, INTENT(IN) :: M
43     ! Septadiagonal matrix A_m
44           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
45     ! Vector b_m
46           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
47     ! Diffusion time-step
48           DOUBLE PRECISION, INTENT(IN) :: lDT
49     !-----------------------------------------------
50     ! Local variables
51     !-----------------------------------------------
52     ! Indices
53           INTEGER :: IJK
54     
55           DOUBLE PRECISION :: APO
56           DOUBLE PRECISION :: lOoDT
57     !-----------------------------------------------
58     
59     
60           lOoDT = 1.0d0/lDT
61     
62           DO IJK = IJKSTART3, IJKEND3
63     
64              IF (FLUID_AT(IJK)) THEN
65     
66                 APO = VOL(IJK)*lOoDT
67     
68     ! Collect the terms
69                 A_M(IJK,0,M) = -(A_M(IJK,E,M)+A_M(IJK,W,M)+&
70                                  A_M(IJK,N,M)+A_M(IJK,S,M)+&
71                                  A_M(IJK,T,M)+A_M(IJK,B,M)+ APO)
72     
73                 S_C = APO*PHI(IJK)
74     
75                 IF(B_M(IJK,M) < S_C .OR. PHI(IJK) == ZERO) THEN
76                    B_M(IJK,M) = B_M(IJK,M) - S_C
77                 ELSE
78                    A_M(IJK,0,M) = A_M(IJK,0,M) - B_M(IJK,M)/PHI(IJK)
79                    B_M(IJK,M) = -S_C
80                 ENDIF
81     
82              ENDIF
83           ENDDO
84     
85     
86           RETURN
87           END SUBROUTINE DIF_PHI_SOURCE_DES
88