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