/[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 1919 - (show annotations) (download) (as text)
Fri Oct 3 05:52:22 2008 UTC (13 years, 8 months ago) by jpye
File MIME type: text/x-csrc
File size: 7230 byte(s)
Implemented (dh/drho)_T
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){
10
11 double rho, T, p, u, h, a, s, cp0;
12
13 double maxerr = 0;
14 double se = 0, sse = 0;
15
16 unsigned i;
17 const unsigned n = ntd;
18
19 fprintf(stderr,"Running through %d test points...\n",n);
20
21 /* Checking CP0 values */
22
23 #define CP0_TEMP(T,RHO,DATA) helmholtz_cp0(T,DATA)
24 fprintf(stderr,"CP0 TESTS\n");
25 for(i=0; i<n;++i){
26 cp0 = td[i].cp0*1e3;
27 ASSERT_TOL(CP0_TEMP, td[i].T+273.15, td[i].rho, d, cp0, cp0*1e-6);
28 }
29 #undef CP0_TEMP
30
31 /* Checking pressure values (proves phir_delta) */
32 fprintf(stderr,"PRESSURE TESTS\n");
33 for(i=16; i<n;++i){
34 T = td[i].T+273.15;
35 rho = td[i].rho;
36 p = td[i].p*1e6;
37 ASSERT_TOL(helmholtz_p, T, rho, d, p, p*1e-3);
38 }
39
40 /* Checking internal energy values (proves phi0_tau, phir_tau) */
41
42 fprintf(stderr,"INTERNAL ENERGY TESTS\n");
43 for(i=0; i<n;++i){
44 T = td[i].T+273.15;
45 rho = td[i].rho;
46 u = td[i].u*1e3;
47 ASSERT_TOL(helmholtz_u, T, rho, d, u, 1e3*u);
48 }
49
50 /* Checking entropy values */
51
52 fprintf(stderr,"ENTROPY TESTS\n");
53 se = 0, sse = 0;
54 for(i=0; i<n;++i){
55 T = td[i].T+273.15;
56 rho = td[i].rho;
57 s = td[i].s*1e3;
58 ASSERT_TOL(helmholtz_s, T, rho, d, s, 1e-6*s);
59 }
60
61 /* checking enthalpy values */
62 fprintf(stderr,"ENTHALPY TESTS\n");
63 se = 0, sse = 0;
64 for(i=0; i<n;++i){
65 T = td[i].T+273.15;
66 rho = td[i].rho;
67 h = td[i].h*1e3;
68 #if 0
69 double err = h - helmholtz_h(T,rho,d);
70 se += err;
71 sse += err*err;
72 fprintf(stderr,"%.12e\t%.12e\t%.12e\t%.12e\n",T,rho,u,err);
73 #else
74 ASSERT_TOL(helmholtz_h, td[i].T+273.15, td[i].rho, d, h, 1E3);
75 #endif
76 }
77 #if 0
78 fprintf(stderr,"average error = %.10e\n",se/n);
79 fprintf(stderr,"sse - n se^2 = %.3e\n",sse - n*se*se);
80 exit(1);
81 #endif
82
83 /* Checking helmholtz energy values */
84
85 fprintf(stderr,"HELMHOLTZ ENERGY TESTS\n");
86 for(i=0; i<n;++i){
87 T = td[i].T+273.15;
88 rho = td[i].rho;
89 a = td[i].a*1e3;
90 //fprintf(stderr,"%.20e\t%.20e\t%.20e\n",T,rho,(a - helmholtz_a(T,rho,d)));
91 ASSERT_TOL(helmholtz_a, T, rho, d, a, a*1e-6);
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 int helm_check_dpdT_rho(const HelmholtzData *d, unsigned ntd, const TestData *td){
141 unsigned i;
142 double T,rho,p,T1,p1,dpdT,dpdT_est,err,se = 0,sse = 0;
143 unsigned n = ntd;
144 double tol = 1e-1;
145
146 double dT = 0.0001 /* finite difference in temperature, in K */;
147
148 fprintf(stderr,"(dp/dT)_rho RESULTS\n\n");
149 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");
150 for(i=0; i<n;++i){
151 T = td[i].T + 273.15;
152 rho = td[i].rho;
153 p = helmholtz_p(T,rho,d);
154 dpdT = helmholtz_dpdT_rho(T,rho,d);
155 assert(!isinf(dpdT));
156 T1 = T + dT;
157 p1 = helmholtz_p(T1, rho, d);
158 dpdT_est = (p1 - p)/dT;
159 err = (dpdT_est - dpdT);
160 se += err;
161 sse += err*err;
162 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 );
163 }
164 fprintf(stderr,"average error = %.10e\n",se/n);
165 fprintf(stderr,"sse - n se^2 = %.3e\n",sse - n*se*se);
166 }
167
168 int helm_check_dpdrho_T(const HelmholtzData *d, unsigned ntd, const TestData *td){
169 unsigned i;
170 double T,rho,p,rho1,p1,dpdrho,dpdrho_est,err,se = 0,sse = 0;
171 unsigned n = ntd;
172 double tol = 1e-1;
173
174 double drho = 0.0001 /* finite difference in temperature, in K */;
175
176 fprintf(stderr,"(dp/drho)_T RESULTS\n\n");
177 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");
178 for(i=0; i<n;++i){
179 T = td[i].T + 273.15;
180 rho = td[i].rho;
181 p = helmholtz_p(T,rho,d);
182 dpdrho = helmholtz_dpdrho_T(T,rho,d);
183 assert(!isinf(dpdrho));
184 rho1 = rho + drho;
185 p1 = helmholtz_p(T, rho1, d);
186 dpdrho_est = (p1 - p)/drho;
187 err = (dpdrho_est - dpdrho);
188 se += err;
189 sse += err*err;
190 fprintf(stderr,"%.12e\t%.12e\t%.12e\t%.12e\t%.12e\t%12.4e\t%12.2f\n",T,rho,p,dpdrho,dpdrho_est,err,err/dpdrho*100);
191 }
192 fprintf(stderr,"average error = %.10e\n",se/n);
193 fprintf(stderr,"sse - n se^2 = %.3e\n",sse - n*se*se);
194 }
195
196
197 int helm_check_dhdT_rho(const HelmholtzData *d, unsigned ntd, const TestData *td){
198 unsigned i;
199 double T,rho,h,T1,h1,dhdT,dhdT_est,err,se = 0,sse = 0;
200 unsigned n = ntd;
201 double tol = 1e-1;
202
203 double dT = 1e-6 /* finite difference in temperature, in K */;
204
205 fprintf(stderr,"(dh/dT)_rho RESULTS\n\n");
206 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");
207 for(i=0; i<n;++i){
208 T = td[i].T + 273.15;
209 rho = td[i].rho;
210 h = helmholtz_h(T,rho,d);
211 dhdT = helmholtz_dhdT_rho(T,rho,d);
212 assert(!isinf(dhdT));
213 T1 = T + dT;
214 h1 = helmholtz_h(T1, rho, d);
215 dhdT_est = (h1 - h)/dT;
216 err = (dhdT_est - dhdT);
217 se += err;
218 sse += err*err;
219 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 );
220 }
221 fprintf(stderr,"average error = %.10e\n",se/n);
222 fprintf(stderr,"sse - n se^2 = %.3e\n",sse - n*se*se);
223 }
224
225
226 int helm_check_dhdrho_T(const HelmholtzData *d, unsigned ntd, const TestData *td){
227 unsigned i;
228 double T,rho,h,rho1,h1,dhdrho,dhdrho_est,err,se = 0,sse = 0;
229 unsigned n = ntd;
230 double tol = 1e-1;
231
232 double drho = 1e-6 /* finite difference in temperature, in K */;
233
234 fprintf(stderr,"(dh/drho)_T RESULTS\n\n");
235 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");
236 for(i=0; i<n;++i){
237 T = td[i].T + 273.15;
238 rho = td[i].rho;
239 h = helmholtz_h(T,rho,d);
240 dhdrho = helmholtz_dhdrho_T(T,rho,d);
241 assert(!isinf(dhdrho));
242 rho1 = rho + drho;
243 h1 = helmholtz_h(T, rho1, d);
244 dhdrho_est = (h1 - h)/drho;
245 err = (dhdrho_est - dhdrho);
246 se += err;
247 sse += err*err;
248 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 );
249 }
250 fprintf(stderr,"average error = %.10e\n",se/n);
251 fprintf(stderr,"sse - n se^2 = %.3e\n",sse - n*se*se);
252 }
253
254

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