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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2682 - (hide 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 jpye 2662
2     #include "../fluids.h"
3     #include "../fprops.h"
4     #include "../solve_ph.h"
5 jpye 2664 #include "../refstate.h"
6 jpye 2682 #include "../sat.h"
7 jpye 2664
8 jpye 2662 #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 jpye 2664 #define TOL_T 1e-8
28 jpye 2662 #define TOL_RHO 1e-3
29    
30    
31 jpye 2682 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 jpye 2662 int main(void){
149 jpye 2664 const PureFluid *P;
150 jpye 2682 //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 jpye 2680 FpropsError err = FPROPS_NO_ERROR;
174 jpye 2662 FluidState S;
175     double T0, rho0, p, h, T, rho;
176    
177 jpye 2682
178    
179 jpye 2662 #define TEST_PH(T1,RHO1) \
180 jpye 2664 err = FPROPS_NO_ERROR;\
181     /*MSG("TEST_PH(T1=%f, RHO1=%f)",T1,RHO1);*/ \
182 jpye 2662 T0 = T1; \
183     rho0 = RHO1; \
184 jpye 2664 /*MSG("setting T,rho");*/\
185 jpye 2662 S = fprops_set_Trho(T0,rho0,P,&err); \
186     assert(!err); \
187 jpye 2664 /*MSG("calc p");*/\
188 jpye 2662 p = fprops_p(S,&err); \
189     assert(!err); \
190 jpye 2664 /*MSG("calc h");*/\
191 jpye 2662 h = fprops_h(S,&err); \
192     assert(!err); \
193 jpye 2664 /*MSG("solving ph");*/\
194 jpye 2662 fprops_solve_ph(p,h,&T,&rho,0,P,&err); \
195     assert(!err); \
196 jpye 2664 MSG("T err: %f",(T-T0));\
197     MSG("rho err: %f",(rho-rho0));\
198 jpye 2662 assert(fabs(T - T0) < TOL_T); \
199     assert(fabs(rho - rho0) < TOL_RHO);
200    
201 jpye 2682 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 jpye 2680 P = fprops_fluid("isohexane","pengrob","J. Chem. Eng. Data, 51");
213     assert(P);
214 jpye 2668
215 jpye 2680 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 jpye 2668 P = fprops_fluid("isohexane","helmholtz","J. Chem. Eng. Data, 51");
223     assert(P);
224 jpye 2680
225     TEST_PH(119.6, 807.530551164909);
226 jpye 2668
227 jpye 2664 P = fprops_fluid("water","helmholtz",NULL);
228     assert(P);
229    
230     TEST_PH(314.4054538268115, 999.7897902747587);
231    
232 jpye 2662 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 jpye 2682 #endif
242 jpye 2662
243 jpye 2680 /* all done? report success */
244 jpye 2662 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