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