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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: LEQ_SOR                                                 C
4     !  Purpose: Solve system of linear system using SOR method             C
5     !           Successive over-relaxation                                 C
6     !                                                                      C
7     !  Author: M. Syamlal                                 Date: 19-AUG-96  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 LEQ_SOR(VNAME, VNO, VAR, A_M, B_M, ITMAX, IER)
18     
19     !-----------------------------------------------
20     ! Modules
21     !-----------------------------------------------
22           USE param
23           USE param1
24           USE matrix
25           USE geometry
26           USE indices
27           USE compar
28           USE sendrecv
29           USE leqsol
30           USE functions
31           IMPLICIT NONE
32     !-----------------------------------------------
33     ! Dummy arguments
34     !-----------------------------------------------
35     ! variable name
36           CHARACTER(LEN=*), INTENT(IN) :: Vname
37     ! variable number (not really used here; see calling subroutine)
38           INTEGER, INTENT(IN) :: VNO
39     ! variable
40     !     e.g., pp_g, epp, rop_g, rop_s, u_g, u_s, v_g, v_s, w_g,
41     !     w_s, T_g, T_s, x_g, x_s, Theta_m, scalar, K_Turb_G,
42     !     e_Turb_G
43           DOUBLE PRECISION, INTENT(INOUT) :: Var(DIMENSION_3)
44     ! Septadiagonal matrix A_m
45           DOUBLE PRECISION, INTENT(INOUT) :: &
46                             A_m(DIMENSION_3, -3:3)
47     ! Vector b_m
48           DOUBLE PRECISION, INTENT(INOUT) :: &
49                             B_m(DIMENSION_3)
50     ! maximum number of iterations (generally leq_it)
51           INTEGER, INTENT(IN) :: ITMAX
52     ! error indicator
53           INTEGER, INTENT(INOUT) :: IER
54     !-----------------------------------------------
55     ! Local parameters
56     !-----------------------------------------------
57     ! OVERRELAXATION FACTOR
58           DOUBLE PRECISION, PARAMETER :: OMEGA = 1.0  !1.2
59           integer :: iidebug
60           parameter( iidebug = 0 )
61     !-----------------------------------------------
62     ! Local Variables
63     !-----------------------------------------------
64     ! Variable
65           DOUBLE PRECISION :: Var_tmp(DIMENSION_3)
66     ! Indices
67           INTEGER :: IJK
68           INTEGER :: ITER
69     
70           DOUBLE PRECISION oAm
71     !-----------------------------------------------
72     
73     !!$omp parallel do private(IJK,OAM)
74           DO IJK = ijkstart3, ijkend3
75              IF(.NOT.IS_ON_myPE_owns(I_OF(IJK),J_OF(IJK), K_OF(IJK))) CYCLE
76     
77              OAM = ONE/A_M(IJK,0)
78              A_M(IJK,0) = ONE
79              A_M(IJK,-2) = A_M(IJK,-2)*OAM
80              A_M(IJK,-1) = A_M(IJK,-1)*OAM
81              A_M(IJK,1) = A_M(IJK,1)*OAM
82              A_M(IJK,2) = A_M(IJK,2)*OAM
83              A_M(IJK,-3) = A_M(IJK,-3)*OAM
84              A_M(IJK,3) = A_M(IJK,3)*OAM
85              B_M(IJK) = B_M(IJK)*OAM
86           ENDDO
87     
88           DO ITER = 1, ITMAX
89              IF (DO_K) THEN
90     
91     !!$omp parallel do private(IJK)
92                 DO IJK = ijkstart3, ijkend3
93                   IF(.NOT.IS_ON_myPE_owns(I_OF(IJK),J_OF(IJK), K_OF(IJK))) CYCLE
94                   VAR_tmp(IJK) = VAR(IJK) + OMEGA*(B_M(IJK)-&
95                      A_M(IJK,-1)*VAR(IM_OF(IJK))-A_M(IJK,1)*VAR(IP_OF(IJK))-&
96                      A_M(IJK,-2)*VAR(JM_OF(IJK))-A_M(IJK,2)*VAR(JP_OF(IJK))-&
97                      A_M(IJK,-3)*VAR(KM_OF(IJK))-A_M(IJK,3)*VAR(KP_OF(IJK))-&
98                      VAR(IJK))
99                 ENDDO
100              ELSE
101     
102     !!$omp parallel do private(IJK)
103                DO IJK = ijkstart3, ijkend3
104                  IF(.NOT.IS_ON_myPE_owns(I_OF(IJK),J_OF(IJK), K_OF(IJK))) CYCLE
105                  VAR_tmp(IJK) = VAR(IJK) + OMEGA*(B_M(IJK)-&
106                       A_M(IJK,-2)*VAR(JM_OF(IJK))-A_M(IJK,2)*VAR(JP_OF(IJK))-&
107                       A_M(IJK,-1)*VAR(IM_OF(IJK))-A_M(IJK,1)*VAR(IP_OF(IJK))-&
108                       VAR(IJK))
109                ENDDO
110              ENDIF
111     
112           call send_recv(var,2)
113           ENDDO
114     
115     !!$omp parallel do private(IJK)
116           DO IJK = ijkstart3, ijkend3
117             VAR(IJK) = VAR_tmp(IJK)
118           ENDDO
119     
120           ITER_TOT(VNO) = ITER
121     
122     
123           RETURN
124           END SUBROUTINE LEQ_SOR
125