File: /nfs/home/0/users/jenkins/mfix.git/model/test_lin_eq.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22 SUBROUTINE TEST_LIN_EQ(A_M, LEQIT, LEQMETHOD, LEQSWEEP, LEQTOL, LEQPC, TEST, IER)
23
24
25
26
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
38
39 INTEGER TEST
40 INTEGER IER
41 DOUBLE PRECISION, DIMENSION(DIMENSION_3,-3:3) :: A_M
42
43
44
45
46
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
56 INTEGER LEQMETHOD, LEQIT
57 CHARACTER(LEN=4) :: LEQSWEEP
58 DOUBLE PRECISION LEQTOL
59 CHARACTER(LEN=4) :: LEQPC
60 REAL :: Harvest
61
62
63
64 CALL RANDOM_SEED
65
66
67
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
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
123
124 CALL SOLVE_LIN_EQ ('Test', 1, X_SOL, Am, Bm, 0, LEQIT, LEQMETHOD, LEQSWEEP, LEQTOL, LEQPC, IER)
125
126
127
128 = 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
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