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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2682 - (show annotations) (download) (as text)
Thu Jan 31 00:41:15 2013 UTC (9 years, 6 months ago) by jpye
File MIME type: text/x-csrc
File size: 6522 byte(s)
Working on better/faster (p,h) test suite.
1
2 #include "../fluids.h"
3 #include "../fprops.h"
4 #include "../solve_ph.h"
5 #include "../refstate.h"
6 #include "../sat.h"
7
8 #include <assert.h>
9 #include <math.h>
10
11 #include "../color.h"
12
13 #define MSG(FMT, ...) \
14 color_on(stderr,ASC_FG_BRIGHTRED);\
15 fprintf(stderr,"%s:%d: ",__FILE__,__LINE__);\
16 color_on(stderr,ASC_FG_BRIGHTBLUE);\
17 fprintf(stderr,"%s: ",__func__);\
18 color_off(stderr);\
19 fprintf(stderr,FMT "\n",##__VA_ARGS__)
20
21 #define ERRMSG(STR,...) \
22 color_on(stderr,ASC_FG_BRIGHTRED);\
23 fprintf(stderr,"ERROR:");\
24 color_off(stderr);\
25 fprintf(stderr," %s:%d:" STR "\n", __func__, __LINE__ ,##__VA_ARGS__)
26
27 #define TOL_T 1e-8
28 #define TOL_RHO 1e-3
29
30
31 void test_ph_array(const PureFluid *P, int nT, int nv, double Trmin, double Trmax, double prskip, int *testedpoints, int *nerr, FpropsError *err){
32 *err = FPROPS_NO_ERROR;
33 *testedpoints = 0;
34 int i,j;
35 double Tmin, Tmax, dT;
36 double vmin, vmax, dlv;
37
38 if(!P){
39 ERRMSG("Pure fluid is not defined");
40 *err = FPROPS_INVALID_REQUEST;
41 return;
42 }
43
44 double Tc = P->data->T_c;
45 double vc = 1./P->data->rho_c;
46 double pc = P->data->p_c;
47 Tmax = Trmax * Tc;
48 Tmin = P->data->T_t;
49 if(0 == Tmin){
50 double psat,rhof, rhog;
51 Tmin = Trmin * Tc;
52 fprops_sat_T(Tmin, &psat, &rhof, &rhog, P, err);
53 if(*err){
54 ERRMSG("Saturation at T = %f failed for '%s",Tmin,P->name);
55 return;
56 }
57 vmin = 1./rhof;
58 vmax = 1./rhog;
59 }else{
60 double pt, rhoft, rhogt;
61 fprops_triple_point(&pt, &rhoft, &rhogt, P, err);
62 if(*err){
63 ERRMSG("Triple point properties failed for '%s",P->name);
64 return;
65 }
66 vmin = 1./rhoft;
67 vmax = 1.2/rhogt;
68 if(vmax > 1e4*vc)vmax=1e4*vc;
69 }
70 /* embarrassingly parallel loops coming up... */
71
72 //#define TOLT 1e-5
73 //#define TOLV 1e-5
74 #define TOLT 1e-1
75 #define TOLV 1e-1
76
77 *nerr = 0;
78
79 FILE *conff = fopen("pherr.ini","w");
80 fprintf(conff,"[main]\n");
81 fprintf(conff,"fluid: %s\n",P->name);
82 fprintf(conff,"type: %s\n",fprops_corr_type(P->type));
83 fprintf(conff,"source: %s\n",P->source);
84 fprintf(conff,"tmin: %e\n",Tmin);
85 fprintf(conff,"tmax: %e\n",Tmax);
86 fprintf(conff,"vmin: %e\n",vmin);
87 fprintf(conff,"vmax: %e\n",vmax);
88 // lower limit for saturation curves, might not be same as tmin in all cases
89 fprintf(conff,"tt: %e\n",Tmin);
90 #define PHERR_DAT "pherr.dat"
91 fprintf(conff,"data: %s\n",PHERR_DAT);
92
93 fclose(conff);
94
95 FILE *plotf = fopen(PHERR_DAT,"w");
96
97 // linear spacing on T, log spacing on v
98 dT = (Tmax - Tmin) / (nT - 1);
99 dlv = (log(vmax) - log(vmin)) / (nv - 1);
100 double T = Tmin, lv = log(vmin);
101 for(i=0; i<nT; ++i, T += dT){
102 fprintf(stderr,">>> T = %6.3f %-40s \r",T,P->name);
103 lv = log(vmin);
104 for(j=0; j<nv; ++j, lv += dlv){
105 double v = exp(lv);
106 double rho = 1./v;
107 *err = FPROPS_NO_ERROR;
108 FluidState S = fprops_set_Trho(T,rho,P,err);
109 if(*err){ERRMSG("Can't set (T,rho)");return;}
110 double p = fprops_p(S,err);
111 if(*err){ERRMSG("Couldn't calculate p");return;}
112
113 if(p / pc > prskip)continue;
114
115 double h = fprops_h(S,err);
116 if(*err){ERRMSG("Couldn't calculate h");return;}
117
118 double T1, rho1;
119 fprops_solve_ph(p,h,&T1,&rho1,0,P,err);
120 (*testedpoints)++;
121 if(*err){
122 /*ERRMSG("Couldn't solve (p,h) at T = %f, rho = %f",T,1./v);*/
123 fprintf(plotf,"%e\t%e\n",T,v);
124 (*nerr)++;
125 continue;
126 }
127
128 double Terr = fabs(T1 - T)/T;
129 double verr = fabs(1./rho1 - v)/v;
130 if(Terr > TOLT || verr > TOLV){
131 //ERRMSG("T or v fail tol at T = %f, rho = %f (Terr %f%%, verr %f%%)",T,rho,Terr*100,verr*100);
132 fprintf(plotf,"%e\t%e\n",T,v);
133 (*nerr)++;
134 continue;
135 }
136 }
137 }
138 fclose(plotf);
139 if(*nerr){
140 ERRMSG("Got %d errors for %s off %d tested points",*nerr,P->name,*testedpoints);
141 *err = FPROPS_NUMERIC_ERROR;
142 }else{
143 ERRMSG("Tested %d points for %s with no failures (T tol %.0e, v tol %.0e)",*testedpoints,P->name,TOLT,TOLV);
144 }
145 }
146
147
148 int main(void){
149 const PureFluid *P;
150 //const char *fluids[] = {"water","toluene","ethanol","isohexane","n_eicosane", NULL};
151 const char *fluids[] = {"isohexane",NULL};
152 const char **fi = fluids;
153 int nfluiderrors = 0;
154 while(*fi){
155 FpropsError err = FPROPS_NO_ERROR;
156 P = fprops_fluid(*fi,"pengrob",NULL);
157 assert(P);
158 int testedpoints;
159 int nerr;
160 test_ph_array(P, /*nT*/300, /*nv*/300, /*Trmin*/0.4, /*Trmax*/1.8, /*prskip*/1.5, &testedpoints, &nerr, &err);
161 if(err)nfluiderrors++;
162 ++fi;
163 }
164
165 if(nfluiderrors){
166 ERRMSG("There were %d fluids with (p,h) errors",nfluiderrors);
167 ERRMSG("Run 'python python/view_ph_results.py' to view results from '%s'",*(fi-1));
168 return 1;
169 }
170
171
172 #if 0
173 FpropsError err = FPROPS_NO_ERROR;
174 FluidState S;
175 double T0, rho0, p, h, T, rho;
176
177
178
179 #define TEST_PH(T1,RHO1) \
180 err = FPROPS_NO_ERROR;\
181 /*MSG("TEST_PH(T1=%f, RHO1=%f)",T1,RHO1);*/ \
182 T0 = T1; \
183 rho0 = RHO1; \
184 /*MSG("setting T,rho");*/\
185 S = fprops_set_Trho(T0,rho0,P,&err); \
186 assert(!err); \
187 /*MSG("calc p");*/\
188 p = fprops_p(S,&err); \
189 assert(!err); \
190 /*MSG("calc h");*/\
191 h = fprops_h(S,&err); \
192 assert(!err); \
193 /*MSG("solving ph");*/\
194 fprops_solve_ph(p,h,&T,&rho,0,P,&err); \
195 assert(!err); \
196 MSG("T err: %f",(T-T0));\
197 MSG("rho err: %f",(rho-rho0));\
198 assert(fabs(T - T0) < TOL_T); \
199 assert(fabs(rho - rho0) < TOL_RHO);
200
201 P = fprops_fluid("water","pengrob","IAPWS");
202 ReferenceState R = {FPROPS_REF_TPF};
203 fprops_set_reference_state(P,&R);
204 MSG("Testing '%s' (type '%d', source '%s')", P->name, P->type, P->source);
205 assert(P);
206 err = FPROPS_NO_ERROR;
207
208 TEST_PH(2.731600000000e+02, 8.618514819566e+02);
209 MSG("p = %f",p);
210 MSG("h = %f",h);
211
212 P = fprops_fluid("isohexane","pengrob","J. Chem. Eng. Data, 51");
213 assert(P);
214
215 TEST_PH(128.4465, 1.);
216 TEST_PH(1.284465e+02, 1/1.693087e-03);
217 TEST_PH(1.284465e+02, 1/2.232836e-03);
218 TEST_PH(1.284465e+02, 1/2.944654e-03);
219 TEST_PH(1.284465e+02, 1/3.883396e-03);
220 TEST_PH(1.284465e+02, 1/5.121405e-03);
221
222 P = fprops_fluid("isohexane","helmholtz","J. Chem. Eng. Data, 51");
223 assert(P);
224
225 TEST_PH(119.6, 807.530551164909);
226
227 P = fprops_fluid("water","helmholtz",NULL);
228 assert(P);
229
230 TEST_PH(314.4054538268115, 999.7897902747587);
231
232 TEST_PH(4.278618181818e+02, 3.591421624513e-03);
233 TEST_PH(3.453541818182e+02, 6.899880592960e-03);
234 TEST_PH(7.166385454545e+02, 6.899880592960e-03);
235 TEST_PH(7.372654545455e+02, 4.092431601778e-03);
236
237 TEST_PH(304.10372142680086, 999.7863801065582);
238 TEST_PH(283.4886584572104, 999.7900620473787);
239 TEST_PH(293.8028752316878, 999.7858245521049);
240
241 #endif
242
243 /* all done? report success */
244 fprintf(stderr,"\n");
245 color_on(stderr,ASC_FG_BRIGHTGREEN);
246 fprintf(stderr,"SUCCESS (%s)",__FILE__);
247 color_off(stderr);
248 fprintf(stderr,"\n");
249 return 0;
250 }
251

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