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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Module name: TEST_LIN_EQ(A_m, LEQIT, LEQMETHOD, LEQSWEEP, LEQTOL, LEQPC, TEST, IER)
4     !  Purpose: Routine for testing the accuracy of linear equation solver C
5     !                                                                      C
6     !                                                                      C
7     !  Author: M. Syamlal                                 Date: 4-JUN-96   C
8     !  Reviewer:                                          Date:            C
9     !                                                                      C
10     !  **** See Solve_energy_eq.f for an example  ****                     C
11     !
12     !  Literature/Document References:                                     C
13     !                                                                      C
14     !  Variables referenced:                                               C
15     !  Variables modified:                                                 C
16     !                                                                      C
17     !  Local variables:                                                    C
18     !                                                                      C
19     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
20     !
21     !
22           SUBROUTINE TEST_LIN_EQ(A_M, LEQIT, LEQMETHOD, LEQSWEEP, LEQTOL, LEQPC,  TEST, IER)
23     !...Translated by Pacific-Sierra Research VAST-90 2.06G5  12:17:31  12/09/98
24     !...Switches: -xf
25     !-----------------------------------------------
26     !   M o d u l e s
27     !-----------------------------------------------
28           USE param
29           USE param1
30           USE matrix
31           USE geometry
32           USE indices
33           USE compar
34           USE functions
35           IMPLICIT NONE
36     !-----------------------------------------------
37     !   D u m m y   A r g u m e n t s
38     !-----------------------------------------------
39           INTEGER TEST  !=0 use the passed A_m; =1 construct a random A_m
40           INTEGER IER
41           DOUBLE PRECISION, DIMENSION(DIMENSION_3,-3:3) :: A_M
42     !-----------------------------------------------
43     !   L o c a l   P a r a m e t e r s
44     !-----------------------------------------------
45     !-----------------------------------------------
46     !   L o c a l   V a r i a b l e s
47     !-----------------------------------------------
48           INTEGER :: IJKERR
49           INTEGER IJK, IpJK, ImJK, IJpK, IJmK, IJKp, IJKm
50           DOUBLE PRECISION, DIMENSION(DIMENSION_3,-3:3) :: Am
51           DOUBLE PRECISION, DIMENSION(DIMENSION_3) :: Bm, X_ACT, X_SOL
52           DOUBLE PRECISION :: ERR, ERRMAX, ERRSUM, XSUM
53           CHARACTER(LEN=80), DIMENSION(8) :: LINE
54     !
55     !                      linear equation solver method and iterations
56           INTEGER          LEQMETHOD, LEQIT
57           CHARACTER(LEN=4) ::   LEQSWEEP
58           DOUBLE PRECISION LEQTOL
59           CHARACTER(LEN=4) ::   LEQPC
60           REAL  :: Harvest
61     
62     !  Initialize the random number generator
63     !
64           CALL RANDOM_SEED
65     !
66     !  Fill the A and x arrays with random numbers, but ensuring that
67     !  the matrix is diagonally dominant
68     !
69           DO IJK = IJKSTART3, IJKEND3
70              CALL RANDOM_NUMBER(HARVEST)
71              X_ACT(IJK) = DBLE(HARVEST) + 1.E-5
72              X_SOL(IJK) = 0.0
73              IF (TEST == 0) THEN
74                 Am(IJK,-3) = A_M(IJK,-3)
75                 Am(IJK,-2) = A_M(IJK,-2)
76                 Am(IJK,-1) = A_M(IJK,-1)
77                 Am(IJK,0) = A_M(IJK,0)
78                 Am(IJK,1) = A_M(IJK,1)
79                 Am(IJK,2) = A_M(IJK,2)
80                 Am(IJK,3) = A_M(IJK,3)
81              ELSE
82                 CALL RANDOM_NUMBER(HARVEST)
83                 Am(IJK,-3) = DBLE(HARVEST)
84                 CALL RANDOM_NUMBER(HARVEST)
85                 Am(IJK,-2) = DBLE(HARVEST)
86                 CALL RANDOM_NUMBER(HARVEST)
87                 Am(IJK,-1) = DBLE(HARVEST)
88                 CALL RANDOM_NUMBER(HARVEST)
89                 Am(IJK,0) = -DBLE(MAX(HARVEST,0.1))*70.
90                 CALL RANDOM_NUMBER(HARVEST)
91                 Am(IJK,1) = DBLE(HARVEST)
92                 CALL RANDOM_NUMBER(HARVEST)
93                 Am(IJK,2) = DBLE(HARVEST)
94                 CALL RANDOM_NUMBER(HARVEST)
95                 Am(IJK,3) = DBLE(HARVEST)
96              ENDIF
97           END DO
98     
99     
100     !
101     !!!$omp  parallel do private( IJK, IpJK, ImJK, IJpK, IJmK, IJKp, IJKm)
102           DO IJK = ijkstart3, ijkend3
103     
104                 ImJK = IM_OF(IJK)
105                 IJmK = JM_OF(IJK)
106                 IpJK = IP_OF(IJK)
107                 IJpK = JP_OF(IJK)
108                 Bm(IJK) = Am(IJK,0)*X_ACT(IJK)
109                 IF(I_OF(IJK) > 1) Bm(IJK) = Bm(IJK) + Am(IJK,W)*X_ACT(ImJK)
110                 IF(I_OF(IJK) < IMAX2) Bm(IJK) = Bm(IJK) +Am(IJK,E)*X_ACT(IpJK)
111                 IF(J_OF(IJK) > 1) Bm(IJK) = Bm(IJK) + Am(IJK,S)*X_ACT(IJmK)
112                 IF(J_OF(IJK) < JMAX2) Bm(IJK) = Bm(IJK) + Am(IJK,N)*X_ACT(IJpK)
113                 IF (DO_K) THEN
114                    IJKm = KM_OF(IJK)
115                    IJKp = KP_OF(IJK)
116                    IF(K_OF(IJK) > 1) Bm(IJK) = Bm(IJK) + Am(IJK,B)*X_ACT(IJKm)
117                    IF(K_OF(IJK) < KMAX2) Bm(IJK) = Bm(IJK) + Am(IJK,T)*X_ACT(IJKp)
118                 ENDIF
119           END DO
120     
121     !
122     !  Solve the linear equation
123     !
124           CALL SOLVE_LIN_EQ ('Test', 1, X_SOL, Am, Bm, 0, LEQIT, LEQMETHOD, LEQSWEEP, LEQTOL, LEQPC, IER)
125     !
126     !  Check the solution
127     !
128           ERRSUM = 0.0
129           XSUM = 0.0
130           ERRMAX = 0.0
131           IJKERR = 0
132           DO IJK = ijkstart3, ijkend3
133              IF (X_ACT(IJK) /= 0.0) THEN
134                 ERR = ABS(X_SOL(IJK)-X_ACT(IJK))/X_ACT(IJK)
135              ELSE IF (X_SOL(IJK) == 0.0) THEN
136                 ERR = 0.0
137              ELSE
138                 ERR = 1.0D32
139              ENDIF
140              IF (ERR > ERRMAX) THEN
141                 ERRMAX = ERR
142                 IJKERR = IJK
143              ENDIF
144              ERRSUM = ERRSUM + ABS(X_SOL(IJK)-X_ACT(IJK))
145              XSUM = XSUM + ABS(X_ACT(IJK))
146     !        print *, ijk, i_of(ijk), j_of(ijk), nint(err * 100.)
147           END DO
148           IF (XSUM /= 0.0) THEN
149              ERR = ERRSUM/XSUM
150           ELSE IF (ERRSUM == 0.0) THEN
151              ERR = 0.0
152           ELSE
153              ERR = 1.0D32
154           ENDIF
155     !
156           IF (ERR < LEQTOL) THEN
157              IER = 0
158              LINE(1) = 'Message: Lin equation solution satisfies tolerance.'
159           ELSE
160              IER = 1
161              LINE(1) = 'Error: Lin equation solution does not satisfy tolerance!'
162           ENDIF
163     !
164           WRITE (LINE(2), *) 'Average normalized error = ', ERR
165           WRITE (LINE(3), *) 'Max normalized error = ', ERRMAX
166           WRITE (LINE(4), *) 'Location of max error = ', IJKERR
167           WRITE (LINE(5), *) 'Xa and Xs @ max error = ', X_ACT(IJKERR), X_SOL(IJKERR)
168           WRITE (LINE(6), *) 'Sample values of actual (Xa) and solution (Xs):'
169           WRITE (LINE(7), '(A,G12.5, A, I6, A, G12.5, A, I6, A, G12.5)') 'Xa(1)=', &
170              X_ACT(1), '  Xa(', IJKMAX2/2, ')=', X_ACT(IJKMAX2/2), '  Xa(', IJKMAX2, ')=', &
171              X_ACT(IJKMAX2)
172           WRITE (LINE(8), '(A,G12.5, A, I6, A, G12.5, A, I6, A, G12.5)') 'Xs(1)=', &
173              X_SOL(1), '  Xs(', IJKMAX2/2, ')=', X_SOL(IJKMAX2/2), '  Xs(', IJKMAX2&
174              , ')=', X_SOL(IJKMAX2)
175     !
176           CALL WRITE_ERROR ('TEST_LIN_EQ', LINE, 8)
177     !
178           RETURN
179           END SUBROUTINE TEST_LIN_EQ
180     
181     
182     
183     
184