/[ascend]/trunk/models/johnpye/fprops/test.c
ViewVC logotype

Contents of /trunk/models/johnpye/fprops/test.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2235 - (show annotations) (download) (as text)
Thu Jul 29 08:45:54 2010 UTC (9 years, 11 months ago) by jpye
File MIME type: text/x-csrc
File size: 10051 byte(s)
we'll satisfy ourselves with a finite difference approximation of helm_resid_deldeldel for the moment, and
move on with the saturation curve problem.
1 #include "test.h"
2
3 #include <assert.h>
4 #include <stdlib.h>
5 #include <stdio.h>
6 #include <math.h>
7
8
9 int helm_run_test_cases(const HelmholtzData *d, unsigned ntd, const TestData *td, int temp_unit){
10
11 double rho, T, p, u, h, a, s, cp0, T_adj = 0;
12
13 if(temp_unit=='C')T_adj = 273.15;
14
15
16 double maxerr = 0;
17 double se = 0, sse = 0;
18
19 unsigned i;
20 const unsigned n = ntd;
21
22 fprintf(stderr,"Running through %d test points...\n",n);
23
24 /* Checking CP0 values */
25
26 #define CP0_TEMP(T,RHO,DATA) helmholtz_cp0(T,DATA)
27 fprintf(stderr,"CP0 TESTS\n");
28 for(i=0; i<n;++i){
29 cp0 = td[i].cp0*1e3;
30 ASSERT_TOL(CP0_TEMP, td[i].T+T_adj, td[i].rho, d, cp0, cp0*1e-6);
31 }
32 #undef CP0_TEMP
33
34 /* Checking pressure values (proves phir_delta) */
35 fprintf(stderr,"PRESSURE TESTS\n");
36 for(i=16; i<n;++i){
37 T = td[i].T+T_adj;
38 rho = td[i].rho;
39 p = td[i].p*1e6;
40 ASSERT_TOL(helmholtz_p, T, rho, d, p, p*1e-3);
41 }
42
43 /* checking enthalpy values */
44 fprintf(stderr,"ENTHALPY TESTS\n");
45 se = 0, sse = 0;
46 for(i=0; i<n;++i){
47 T = td[i].T+T_adj;
48 rho = td[i].rho;
49 h = td[i].h*1e3;
50 ASSERT_TOL(helmholtz_h, td[i].T+T_adj, td[i].rho, d, h, 1e-3*h);
51 }
52
53 /* Checking internal energy values (proves phi0_tau, phir_tau) */
54 fprintf(stderr,"INTERNAL ENERGY TESTS\n");
55 for(i=0; i<n;++i){
56 T = td[i].T+T_adj;
57 rho = td[i].rho;
58 u = td[i].u*1e3;
59 ASSERT_TOL(helmholtz_u, T, rho, d, u, 1e3*u);
60 }
61
62 /* Checking entropy values */
63
64 fprintf(stderr,"ENTROPY TESTS\n");
65 se = 0, sse = 0;
66 for(i=0; i<n;++i){
67 T = td[i].T+T_adj;
68 rho = td[i].rho;
69 s = td[i].s*1e3;
70 ASSERT_TOL(helmholtz_s, T, rho, d, s, 1e-5*s);
71 }
72
73
74 /* Checking helmholtz energy values */
75
76 fprintf(stderr,"HELMHOLTZ ENERGY TESTS\n");
77 for(i=0; i<n;++i){
78 T = td[i].T+T_adj;
79 rho = td[i].rho;
80 a = td[i].a*1e3;
81 //fprintf(stderr,"%.20e\t%.20e\t%.20e\n",T,rho,(a - helmholtz_a(T,rho,d)));
82 ASSERT_TOL(helmholtz_a, T, rho, d, a, a*1e-3 );
83 }
84
85 fprintf(stderr,"Tests completed OK (maximum error = %0.5f%%)\n",maxerr);
86 exit(0);
87 }
88
89 int helm_check_u(const HelmholtzData *d, unsigned ntd, const TestData *td){
90 unsigned i;
91 double T,rho,u,err,se = 0,sse = 0;
92 unsigned n = ntd;
93 double tol = 1e-1;
94
95 fprintf(stderr,"INTERNAL ENERGY VALUES\n\n");
96 fprintf(stderr,"%-18s\t%-18s\t%-18s\t%-18s\t%-18s\n","T","rho","u","du","%err");
97 for(i=0; i<n;++i){
98 T = td[i].T+273.15;
99 rho = td[i].rho;
100 u = td[i].u*1e3;
101 err = u - helmholtz_u(T,rho,d);
102 se += err;
103 sse += err*err;
104 fprintf(stderr,"%.12e\t%.12e\t%.12e\t%.12e\t%.6f\n",T,rho,u,err,err/u*100.);
105 }
106 fprintf(stderr,"average error = %.10e\n",se/n);
107 fprintf(stderr,"sse - n se^2 = %.3e\n",sse - n*se*se);
108 }
109
110 int helm_check_s(const HelmholtzData *d, unsigned ntd, const TestData *td){
111 unsigned i;
112 double T,rho,s,err,se = 0,sse = 0;
113 unsigned n = ntd;
114 double tol = 1e-1;
115
116 fprintf(stderr,"ENTROPY VALUES\n\n");
117 fprintf(stderr,"%-18s\t%-18s\t%-18s\t%-18s\t%-18s\n","T","rho","s","ds","%err");
118 for(i=0; i<n;++i){
119 T = td[i].T + 273.15;
120 rho = td[i].rho;
121 s = td[i].s*1e3;
122 err = s - helmholtz_s(T,rho,d);
123 se += err;
124 sse += err*err;
125 fprintf(stderr,"%.12e\t%.12e\t%.12e\t%.12e\t%.6f\n",T,rho,s,err,err/s*100.);
126 }
127 fprintf(stderr,"average error = %.10e\n",se/n);
128 fprintf(stderr,"sse - n se^2 = %.3e\n",sse - n*se*se);
129 }
130
131 /* check derivatives of p */
132
133 int helm_check_dpdT_rho(const HelmholtzData *d, unsigned ntd, const TestData *td){
134 unsigned i;
135 double T,rho,p,T1,p1,dpdT,dpdT_est,err,se = 0,sse = 0;
136 unsigned n = ntd;
137 double tol = 1e-1;
138
139 double dT = 0.0001 /* finite difference in temperature, in K */;
140
141 fprintf(stderr,"(dp/dT)_rho RESULTS\n\n");
142 fprintf(stderr,"%-18s\t%-18s\t%-18s\t%-18s\t%-18s\t%12s\t%12s\n","T","rho","p","dp/dT","dp/dT est","err","%err");
143 for(i=0; i<n;++i){
144 T = td[i].T + 273.15;
145 rho = td[i].rho;
146 p = helmholtz_p(T,rho,d);
147 dpdT = helmholtz_dpdT_rho(T,rho,d);
148 assert(!isinf(dpdT));
149 T1 = T + dT;
150 p1 = helmholtz_p(T1, rho, d);
151 dpdT_est = (p1 - p)/dT;
152 err = (dpdT_est - dpdT);
153 se += err;
154 sse += err*err;
155 fprintf(stderr,"%.12e\t%.12e\t%.12e\t%.12e\t%.12e\t%12.4e\t%12.2e\n",T,rho,p,dpdT,dpdT_est,err,err/dpdT*100 );
156 }
157 fprintf(stderr,"average error = %.10e\n",se/n);
158 fprintf(stderr,"sse - n se^2 = %.3e\n",sse - n*se*se);
159 }
160
161 int helm_check_dpdrho_T(const HelmholtzData *d, unsigned ntd, const TestData *td){
162 unsigned i;
163 double T,rho,p,rho1,p1,dpdrho,dpdrho_est,err,se = 0,sse = 0;
164 unsigned n = ntd;
165 double tol = 1e-1;
166
167 double drho = 0.0001 /* finite difference in temperature, in K */;
168
169 fprintf(stderr,"(dp/drho)_T RESULTS\n\n");
170 fprintf(stderr,"%-18s\t%-18s\t%-18s\t%-18s\t%-18s\t%12s\t%12s\n","T","rho","p","dp/drho","dp/drho est","err","%err");
171 for(i=0; i<n;++i){
172 T = td[i].T + 273.15;
173 rho = td[i].rho;
174 p = helmholtz_p(T,rho,d);
175 dpdrho = helmholtz_dpdrho_T(T,rho,d);
176 assert(!isinf(dpdrho));
177 rho1 = rho + drho;
178 p1 = helmholtz_p(T, rho1, d);
179 dpdrho_est = (p1 - p)/drho;
180 err = (dpdrho_est - dpdrho);
181 se += err;
182 sse += err*err;
183 fprintf(stderr,"%.12e\t%.12e\t%.12e\t%.12e\t%.12e\t%12.4e\t%12.6f\n",T,rho,p,dpdrho,dpdrho_est,err,err/dpdrho*100);
184 }
185 fprintf(stderr,"average error = %.10e\n",se/n);
186 fprintf(stderr,"sse - n se^2 = %.3e\n",sse - n*se*se);
187 }
188
189
190 int helm_check_d2pdrho2_T(const HelmholtzData *d, unsigned ntd, const TestData *td){
191 unsigned i;
192 double T,rho,p,dpdrho,rho1,dpdrho1,d2pdrho2,d2pdrho2_est,err,se = 0,sse = 0;
193 unsigned n = ntd;
194 double tol = 1e-1;
195
196 double drho = 0.0001 /* finite difference in temperature, in K */;
197
198 fprintf(stderr,"\n(d2p/drho2)_T RESULTS\n\n");
199 fprintf(stderr,"%-18s\t%-18s\t%-18s\t%-18s\t%-18s\t%12s\t%12s\t%12s\n","T","rho","p","dp/drho","d2p/drho2 est","d2p/drho2 calc","err","rel err");
200 for(i=0; i<n;++i){
201 T = td[i].T + 273.15;
202 rho = td[i].rho;
203 p = helmholtz_p(T,rho,d);
204 dpdrho = helmholtz_dpdrho_T(T,rho,d);
205 d2pdrho2 = helmholtz_d2pdrho2_T(T, rho, d);
206 assert(!isinf(d2pdrho2));
207 rho1 = rho + drho;
208 dpdrho1 = helmholtz_dpdrho_T(T, rho1, d);
209 d2pdrho2_est = (dpdrho1 - dpdrho)/drho;
210 err = (d2pdrho2_est - d2pdrho2);
211 se += err;
212 sse += err*err;
213 fprintf(stderr,"%.12e\t%.12e\t%.12e\t%.12e\t%.12e\t%12.4e\t%12.2e\t%12.3e\n",T,rho,p,dpdrho,d2pdrho2_est,d2pdrho2,err,err/d2pdrho2_est);
214 }
215 fprintf(stderr,"average error = %.10e\n",se/n);
216 fprintf(stderr,"sse - n se^2 = %.3e\n",sse - n*se*se);
217 }
218
219 /* check derivatives of h */
220
221 int helm_check_dhdT_rho(const HelmholtzData *d, unsigned ntd, const TestData *td){
222 unsigned i;
223 double T,rho,h,T1,h1,dhdT,dhdT_est,err,se = 0,sse = 0;
224 unsigned n = ntd;
225 double tol = 1e-1;
226
227 double dT = 1e-6 /* finite difference in temperature, in K */;
228
229 fprintf(stderr,"(dh/dT)_rho RESULTS\n\n");
230 fprintf(stderr,"%-18s\t%-18s\t%-18s\t%-18s\t%-18s\t%12s\t%12s\n","T","rho","h","dh/dT","dh/dT est","err","%err");
231 for(i=0; i<n;++i){
232 T = td[i].T + 273.15;
233 rho = td[i].rho;
234 h = helmholtz_h(T,rho,d);
235 dhdT = helmholtz_dhdT_rho(T,rho,d);
236 assert(!isinf(dhdT));
237 T1 = T + dT;
238 h1 = helmholtz_h(T1, rho, d);
239 dhdT_est = (h1 - h)/dT;
240 err = (dhdT_est - dhdT);
241 se += err;
242 sse += err*err;
243 fprintf(stderr,"%.12e\t%.12e\t%.12e\t%.12e\t%.12e\t%12.4e\t%12.2e\n",T,rho,h,dhdT,dhdT_est,err,err/dhdT*100 );
244 }
245 fprintf(stderr,"average error = %.10e\n",se/n);
246 fprintf(stderr,"sse - n se^2 = %.3e\n",sse - n*se*se);
247 }
248
249
250 int helm_check_dhdrho_T(const HelmholtzData *d, unsigned ntd, const TestData *td){
251 unsigned i;
252 double T,rho,h,rho1,h1,dhdrho,dhdrho_est,err,se = 0,sse = 0;
253 unsigned n = ntd;
254
255 double drho = 1e-6 /* finite difference in temperature, in K */;
256
257 fprintf(stderr,"(dh/drho)_T RESULTS\n\n");
258 fprintf(stderr,"%-18s\t%-18s\t%-18s\t%-18s\t%-18s\t%12s\t%12s\n","T","rho","h","dh/drho","dh/drho est","err","%err");
259 for(i=0; i<n;++i){
260 T = td[i].T + 273.15;
261 rho = td[i].rho;
262 h = helmholtz_h(T,rho,d);
263 dhdrho = helmholtz_dhdrho_T(T,rho,d);
264 assert(!isinf(dhdrho));
265 rho1 = rho + drho;
266 h1 = helmholtz_h(T, rho1, d);
267 dhdrho_est = (h1 - h)/drho;
268 err = (dhdrho_est - dhdrho);
269 se += err;
270 sse += err*err;
271 fprintf(stderr,"%.12e\t%.12e\t%.12e\t%.12e\t%.12e\t%12.4e\t%12.2e\n",T,rho,h,dhdrho,dhdrho_est,err,err/dhdrho*100 );
272 }
273 fprintf(stderr,"average error = %.10e\n",se/n);
274 fprintf(stderr,"sse - n se^2 = %.3e\n",sse - n*se*se);
275 }
276
277 /* check derivatives of u */
278
279 int helm_check_dudT_rho(const HelmholtzData *d, unsigned ntd, const TestData *td){
280 unsigned i;
281 double T,rho,u,T1,u1,dudT,dudT_est,err,se = 0,sse = 0;
282 unsigned n = ntd;
283
284 double dT = 1e-3 /* finite difference in temperature, in K */;
285
286 fprintf(stderr,"(du/dT)_rho RESULTS\n\n");
287 fprintf(stderr,"%-18s\t%-18s\t%-18s\t%-18s\t%-18s\t%12s\t%12s\n","T","rho","h","du/dT","du/dT est","err","%err");
288 for(i=0; i<n;++i){
289 T = td[i].T + 273.15;
290 rho = td[i].rho;
291 u = helmholtz_u(T,rho,d);
292 dudT = helmholtz_dudT_rho(T,rho,d);
293 assert(!isinf(dudT));
294 T1 = T + dT;
295 u1 = helmholtz_u(T1, rho, d);
296 dudT_est = (u1 - u)/dT;
297 err = (dudT_est - dudT);
298 se += err;
299 sse += err*err;
300 fprintf(stderr,"%.12e\t%.12e\t%.12e\t%.12e\t%.12e\t%12.4e\t%12.2e\n",T,rho,u,dudT,dudT_est,err,err/dudT*100 );
301 }
302 fprintf(stderr,"average error = %.10e\n",se/n);
303 fprintf(stderr,"sse - n se^2 = %.3e\n",sse - n*se*se);
304 }
305
306
307 int helm_check_dudrho_T(const HelmholtzData *d, unsigned ntd, const TestData *td){
308 unsigned i;
309 double T,rho,u,rho1,u1,dudrho,dudrho_est,err,se = 0,sse = 0;
310 unsigned n = ntd;
311 double tol = 1e-1;
312
313 double drho = 1e-6 /* finite difference in temperature, in K */;
314
315 fprintf(stderr,"(du/drho)_T RESULTS\n\n");
316 fprintf(stderr,"%-18s\t%-18s\t%-18s\t%-18s\t%-18s\t%12s\t%12s\n","T","rho","h","du/drho","du/drho est","err","%err");
317 for(i=0; i<n;++i){
318 T = td[i].T + 273.15;
319 rho = td[i].rho;
320 u = helmholtz_u(T,rho,d);
321 dudrho = helmholtz_dudrho_T(T,rho,d);
322 assert(!isinf(dudrho));
323 rho1 = rho + drho;
324 u1 = helmholtz_u(T, rho1, d);
325 dudrho_est = (u1 - u)/drho;
326 err = (dudrho_est - dudrho);
327 se += err;
328 sse += err*err;
329 fprintf(stderr,"%.12e\t%.12e\t%.12e\t%.12e\t%.12e\t%12.4e\t%12.2e\n",T,rho,u,dudrho,dudrho_est,err,err/dudrho*100 );
330 }
331 fprintf(stderr,"average error = %.10e\n",se/n);
332 fprintf(stderr,"sse - n se^2 = %.3e\n",sse - n*se*se);
333 }
334
335
336
337
338
339
340

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22