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

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