1#include <stdio.h>
2#include <math.h>
3#include "applu_share.h"
4#define class _class_
5
6#define xcr(m) xcr[m-1]
7#define xce(m) xce[m-1]
8#define xcrref(m) xcrref[m-1]
9#define xceref(m) xceref[m-1]
10#define xcrdif(m) xcrdif[m-1]
11#define xcedif(m) xcedif[m-1]
12
13void verify(double *xcr, double *xce, double *xcip, char *classp) {
14
15        double xcrref[5],xceref[5],xciref, xcrdif[5], xcedif[5], xcidif,
16               epsilon, dtref;
17        int m;
18
19//c   tolerance level
20
21        double xci = *xcip;
22
23        epsilon = 1.0e-8;
24        char class = 'U';
25        int verified = 1;
26
27        for (m=1; m<=5; m++) {
28           xcrref(m) = 1.0;
29           xceref(m) = 1.0;
30        }
31        xciref = 1.0;
32
33        if ((nx0 == 12) && (ny0 == 12) && (nz0 == 12) && (itmax == 50 )) {
34
35           class = 'S';
36           dtref = 5.0e-1;
37//c---------------------------------------------------------------------
38//c   Reference values of RMS-norms of residual, for the (12X12X12) grid,
39//c   after 50 time steps, with  DT = 5.0d-01
40//---------------------------------------------------------------------
41         xcrref(1) = 1.6196343210976702e-02;
42         xcrref(2) = 2.1976745164821318e-03;
43         xcrref(3) = 1.5179927653399185e-03;
44         xcrref(4) = 1.5029584435994323e-03;
45         xcrref(5) = 3.4264073155896461e-02;
46
47//c---------------------------------------------------------------------
48//c   Reference values of RMS-norms of solution error, for the (12X12X12) grid,
49//c   after 50 time steps, with  DT = 5.0e-01;
50//c---------------------------------------------------------------------
51         xceref(1) = 6.4223319957960924e-04;
52         xceref(2) = 8.4144342047347926e-05;
53         xceref(3) = 5.8588269616485186e-05;
54         xceref(4) = 5.8474222595157350e-05;
55         xceref(5) = 1.3103347914111294e-03;
56
57//c---------------------------------------------------------------------
58//c   Reference value of surface integral, for the (12X12X12) grid,
59//c   after 50 time steps, with DT = 5.0d-01
60//c---------------------------------------------------------------------
61         xciref = 7.8418928865937083;
62
63        } else
64        if ((nx0 == 33) && (ny0 == 33) && (nz0 == 33) && (itmax == 300)) {
65
66           class = 'W';
67           dtref = 1.5e-3;
68//c---------------------------------------------------------------------
69//c   Reference values of RMS-norms of residual, for the (33x33x33) grid,
70//c   after 300 time steps, with  DT = 1.5d-3
71//c---------------------------------------------------------------------
72           xcrref(1) =   0.1236511638192e+02;
73           xcrref(2) =   0.1317228477799e+01;
74           xcrref(3) =   0.2550120713095e+01;
75           xcrref(4) =   0.2326187750252e+01;
76           xcrref(5) =   0.2826799444189e+02;
77
78//c---------------------------------------------------------------------
79//c   Reference values of RMS-norms of solution error, for the (33X33X33) grid,
80//c---------------------------------------------------------------------
81           xceref(1) =   0.4867877144216e+00;
82           xceref(2) =   0.5064652880982e-01;
83           xceref(3) =   0.9281818101960e-01;
84           xceref(4) =   0.8570126542733e-01;
85           xceref(5) =   0.1084277417792e+01;
86
87//c---------------------------------------------------------------------
88//c   Reference value of surface integral, for the (33X33X33) grid,
89//c   after 300 time steps, with  DT = 1.5e-3;
90//c---------------------------------------------------------------------
91           xciref    =   0.1161399311023e+02;
92
93        } else
94        if ((nx0 == 64) && (ny0 == 64) && (nz0 == 64) && (itmax == 250)) {
95
96           class = 'A';
97           dtref = 2.0e+0;
98//c---------------------------------------------------------------------
99//c   Reference values of RMS-norms of residual, for the (64X64X64) grid,
100//c   after 250 time steps, with  DT = 2.0e+00;
101//c---------------------------------------------------------------------
102         xcrref(1) = 7.7902107606689367e+02;
103         xcrref(2) = 6.3402765259692870e+01;
104         xcrref(3) = 1.9499249727292479e+02;
105         xcrref(4) = 1.7845301160418537e+02;
106         xcrref(5) = 1.8384760349464247e+03;
107
108//c---------------------------------------------------------------------
109//c   Reference values of RMS-norms of solution error, for the (64X64X64) grid,
110//c   after 250 time steps, with  DT = 2.0e+00;
111//c---------------------------------------------------------------------
112         xceref(1) = 2.9964085685471943e+01;
113         xceref(2) = 2.8194576365003349e+00;
114         xceref(3) = 7.3473412698774742e+00;
115         xceref(4) = 6.7139225687777051e+00;
116         xceref(5) = 7.0715315688392578e+01;
117
118//c---------------------------------------------------------------------
119//c   Reference value of surface integral, for the (64X64X64) grid,
120//c   after 250 time steps, with DT = 2.0e+00;
121//c---------------------------------------------------------------------
122         xciref = 2.6030925604886277e+01;
123
124        } else
125        if ((nx0 == 102) && (ny0 == 102) && (nz0 == 102) && (itmax == 250)) {
126
127           class = 'B';
128           dtref = 2.0e+0;
129
130//c---------------------------------------------------------------------
131//c   Reference values of RMS-norms of residual, for the (102X102X102) grid,
132//c   after 250 time steps, with  DT = 2.0e+00;
133//c---------------------------------------------------------------------
134         xcrref(1) = 3.5532672969982736e+03;
135         xcrref(2) = 2.6214750795310692e+02;
136         xcrref(3) = 8.8333721850952190e+02;
137         xcrref(4) = 7.7812774739425265e+02;
138         xcrref(5) = 7.3087969592545314e+03;
139
140//c---------------------------------------------------------------------
141//c   Reference values of RMS-norms of solution error, for the (102X102X102)
142//c   grid, after 250 time steps, with  DT = 2.0e+00;
143//c---------------------------------------------------------------------
144         xceref(1) = 1.1401176380212709e+02;
145         xceref(2) = 8.1098963655421574e+00;
146         xceref(3) = 2.8480597317698308e+01;
147         xceref(4) = 2.5905394567832939e+01;
148         xceref(5) = 2.6054907504857413e+02;
149
150//c---------------------------------------------------------------------
151//c   Reference value of surface integral, for the (102X102X102) grid,
152//c   after 250 time steps, with DT = 2.0e+00;
153//c---------------------------------------------------------------------
154         xciref = 4.7887162703308227e+01;
155
156        } else
157        if ((nx0 == 162) && (ny0 == 162) && (nz0 == 162) && (itmax == 250)) {
158
159           class = 'C';
160           dtref = 2.0e+0;
161
162//c---------------------------------------------------------------------
163//c   Reference values of RMS-norms of residual, for the (162X162X162) grid,
164//c   after 250 time steps, with  DT = 2.0e+00;
165//c---------------------------------------------------------------------
166         xcrref(1) = 1.03766980323537846e+04;
167         xcrref(2) = 8.92212458801008552e+02;
168         xcrref(3) = 2.56238814582660871e+03;
169         xcrref(4) = 2.19194343857831427e+03;
170         xcrref(5) = 1.78078057261061185e+04;
171
172//c---------------------------------------------------------------------
173//c   Reference values of RMS-norms of solution error, for the (162X162X162)
174//c   grid, after 250 time steps, with  DT = 2.0e+00;
175//c---------------------------------------------------------------------
176         xceref(1) = 2.15986399716949279e+02;
177         xceref(2) = 1.55789559239863600e+01;
178         xceref(3) = 5.41318863077207766e+01;
179         xceref(4) = 4.82262643154045421e+01;
180         xceref(5) = 4.55902910043250358e+02;
181
182//c---------------------------------------------------------------------
183//c   Reference value of surface integral, for the (162X162X162) grid,
184//c   after 250 time steps, with DT = 2.0e+00;
185//c---------------------------------------------------------------------
186         xciref = 6.66404553572181300e+01;
187
188        } else
189        if ((nx0 == 408) && (ny0 == 408) && (nz0 == 408) && (itmax == 300)) {
190
191           class = 'D';
192           dtref = 1.0e+0;
193
194//c---------------------------------------------------------------------
195//c   Reference values of RMS-norms of residual, for the (408X408X408) grid,
196//c   after 300 time steps, with  DT = 1.0e+00;
197//c---------------------------------------------------------------------
198         xcrref(1) = 0.4868417937025e+05;
199         xcrref(2) = 0.4696371050071e+04;
200         xcrref(3) = 0.1218114549776e+05;
201         xcrref(4) = 0.1033801493461e+05;
202         xcrref(5) = 0.7142398413817e+05;
203
204//c---------------------------------------------------------------------
205//c   Reference values of RMS-norms of solution error, for the (408X408X408)
206//c   grid, after 300 time steps, with  DT = 1.0e+00;
207//c---------------------------------------------------------------------
208         xceref(1) = 0.3752393004482e+03;
209         xceref(2) = 0.3084128893659e+02;
210         xceref(3) = 0.9434276905469e+02;
211         xceref(4) = 0.8230686681928e+02;
212         xceref(5) = 0.7002620636210e+03;
213
214//c---------------------------------------------------------------------
215//c   Reference value of surface integral, for the (408X408X408) grid,
216//c   after 300 time steps, with DT = 1.0e+00;
217//c---------------------------------------------------------------------
218         xciref =    0.8334101392503e+02;
219
220        } else
221        if ((nx0 == 1020) && (ny0 == 1020) && (nz0 == 1020) && (itmax == 300) ) {
222
223           class = 'E';
224           dtref = 0.5e+0;
225
226//c---------------------------------------------------------------------
227//c   Reference values of RMS-norms of residual, for the (1020X1020X1020) grid,
228//c   after 300 time steps, with  DT = 0.5e+00;
229//c---------------------------------------------------------------------
230         xcrref(1) = 0.2099641687874e+06;
231         xcrref(2) = 0.2130403143165e+05;
232         xcrref(3) = 0.5319228789371e+05;
233         xcrref(4) = 0.4509761639833e+05;
234         xcrref(5) = 0.2932360006590e+06;
235
236//c---------------------------------------------------------------------
237//c   Reference values of RMS-norms of solution error, for the (1020X1020X1020)
238//c   grid, after 300 time steps, with  DT = 0.5e+00;
239//c---------------------------------------------------------------------
240         xceref(1) = 0.4800572578333e+03;
241         xceref(2) = 0.4221993400184e+02;
242         xceref(3) = 0.1210851906824e+03;
243         xceref(4) = 0.1047888986770e+03;
244         xceref(5) = 0.8363028257389e+03;
245
246//c---------------------------------------------------------------------
247//c   Reference value of surface integral, for the (1020X1020X1020) grid,
248//c   after 300 time steps, with DT = 0.5e+00;
249//c---------------------------------------------------------------------
250//         xciref =    0.9512163272273e+02;
251
252        } else  verified = 0;
253
254//c---------------------------------------------------------------------
255//c    verification test for residuals if gridsize is one of
256//c    the defined grid sizes above (class != 'U')
257//c---------------------------------------------------------------------
258
259//c---------------------------------------------------------------------
260//c    Compute the difference of solution values and the known reference values.
261//c---------------------------------------------------------------------
262        for (m=1; m<=5; m++) {
263
264           xcrdif(m) = fabs((xcr(m)-xcrref(m))/xcrref(m));
265           xcedif(m) = fabs((xce(m)-xceref(m))/xceref(m));
266
267        }
268        xcidif = fabs((xci - xciref)/xciref);
269
270//c---------------------------------------------------------------------
271//c    Output the comparison of computed results to known cases.
272//c---------------------------------------------------------------------
273
274        if (class != 'U') {
275           printf(" Verification being performed for class %c\n", class);
276           printf(" Accuracy setting for epsilon = %20.13e\n", epsilon);
277           verified = (fabs(dt-dtref) <= epsilon);
278           if (!verified) {
279              class = 'U';
280              printf(" DT does not match the reference value of %15.8e\n",
281                       dtref);
282           }
283        } else printf(" Unknown class\n");
284
285        if (class != 'U') printf(" Comparison of RMS-norms of residual\n");
286        else  printf(" RMS-norms of residual\n");
287
288        for (m=1; m<=5; m++) {
289           if (class == 'U')
290              printf("       %2d  %20.13e\n",  m, xcr(m));
291           else  if (xcrdif(m) <= epsilon) {
292              printf("         %2d  %20.13e  %20.13e  %20.13e\n",
293                      m,xcr(m),xcrref(m),xcrdif(m));
294           } else {
295              verified = 0;
296              printf("FAILURE: %2d  %20.13e  %20.13e  %20.13e\n",
297                      m,xcr(m),xcrref(m),xcrdif(m));
298           }
299        }
300
301        if (class != 'U') printf(" Comparison of RMS-norms of solution error\n");
302         else printf(" RMS-norms of solution error\n");
303
304        for (m=1; m<=5; m++) {
305           if (class == 'U')
306              printf("       %2d  %20.13e\n",  m, xce(m));
307           else  if (xcedif(m) <= epsilon) {
308              printf("         %2d  %20.13e  %20.13e  %20.13e\n",
309                      m,xce(m),xceref(m),xcedif(m));
310           } else {
311              verified = 0;
312              printf("FAILURE: %2d  %20.13e  %20.13e  %20.13e\n",
313                      m,xce(m),xceref(m),xcedif(m));
314           }
315        }
316
317        if (class != 'U') printf(" Comparison of surface integral\n");
318        else  printf(" Surface integral\n");
319
320        if (class == 'U') {
321           printf("             %20.13e\n", xci);
322        } else  if (xcidif <= epsilon) {
323           printf("             %20.13e  %20.13e  %20.13e\n", xci, xciref, xcidif);
324        } else {
325           verified = 0;
326           printf("FAILURE:     %20.13e  %20.13e  %20.13e\n", xci, xciref, xcidif);
327        }
328
329
330        if (class == 'U') {
331           printf(" No reference values provided\n");
332           printf(" No verification performed\n");
333        } else  if (verified) {
334           printf(" Verification Successful\n");
335        } else {
336           printf(" Verification failed\n");
337        }
338
339        *classp    = class;
340
341        return;
342}
343
344