/[ascend]/branches/sid/models/johnpye/fprops/test/test_ttse.c
ViewVC logotype

Contents of /branches/sid/models/johnpye/fprops/test/test_ttse.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3025 - (show annotations) (download) (as text)
Thu Jul 23 18:43:12 2015 UTC (2 years, 11 months ago) by sid
File MIME type: text/x-csrc
File size: 6009 byte(s)


1
2
3 #include "../fluids.h"
4 #include "../fprops.h"
5 #include "../solve_ph.h"
6 #include "../refstate.h"
7 #include "../sat.h"
8 #include "../visc.h"
9 #include "../ttse.h"
10
11 #include <assert.h>
12 #include <math.h>
13 #include <stdio.h>
14 #include <stdlib.h>
15 #include <time.h>
16
17 #include "../color.h"
18
19 #define MSG(FMT, ...) \
20 color_on(stderr,ASC_FG_BRIGHTRED);\
21 fprintf(stderr,"%s:%d: ",__FILE__,__LINE__);\
22 color_on(stderr,ASC_FG_BRIGHTBLUE);\
23 fprintf(stderr,"%s: ",__func__);\
24 color_off(stderr);\
25 fprintf(stderr,FMT "\n",##__VA_ARGS__)
26
27 #define ERRMSG(STR,...) \
28 color_on(stderr,ASC_FG_BRIGHTRED);\
29 fprintf(stderr,"ERROR:");\
30 color_off(stderr);\
31 fprintf(stderr," %s:%d:" STR "\n", __func__, __LINE__ ,##__VA_ARGS__)
32
33 #define TOL_T 1e-8
34 #define TOL_RHO 1e-3
35
36 int main(void){
37 PureFluid *P;
38 FpropsError err;
39 const char *helmfluids[] = { "water"};
40 //const int n = sizeof(helmfluids)/sizeof(char *);
41 int i;
42
43
44 MSG("Which Fluid? --> %s",helmfluids[0]);
45
46 P = (PureFluid *)fprops_fluid(helmfluids[0],"helmholtz",NULL);
47
48
49 MSG("Triple Point Tt %f & Crit Temp %f",P->data->T_t,P->data->T_c);
50
51
52 MSG("Comparing Helmholtz vs TTSE");
53
54
55 double rho = 1000;
56
57 //Plot in Mathematica Saturation dome
58
59 printf("{");
60 for(i=0;i<NSAT;i++)
61 {
62 double Tt = P->data->T_t;
63 double Tc = P->data->T_c;
64 double dt2p = (Tc - Tt)/NSAT;
65 double T = Tt + (i+0.1)*dt2p;
66 double psat, rhof,rhog;
67 fprops_sat_T(T,&psat,&rhof,&rhog,P,&err);
68
69 int j = (int)round( ((T - Tt)/(Tc - Tt))*(NSAT) );
70 double delt = T - ( Tt + j*dt2p);
71 double rhofT,rhogT;
72 rhofT = P->table->satFRho[j] + delt*P->table->satFdRhodt[j]+ 0.5*delt*delt*P->table->satFd2RhodT2[j];
73 rhogT = P->table->satGRho[j] + delt*P->table->satGdRhodt[j]+ 0.5*delt*delt*P->table->satGd2RhodT2[j];
74
75 MSG("%f %f %f %f %f %f %f", delt,rhof,rhofT,rhog,rhogT, 100*(rhof-rhofT)/rhof,100*(rhog-rhogT)/rhog );
76 // if(i!= SPOINTS)
77 // printf("{%f, %f, %f},\n",T, rhof,rhog);
78 // else
79 // printf("{%f, %f, %f}};\n",T, rhof,rhog);
80 }
81
82 exit(1);
83
84
85
86
87
88 /*****************************************Two Phase Table Testing*****************************************/
89 /*
90 double Tt = P->data->T_t;
91 double Tc = P->data->T_c;
92 double t = Tc - 53.341233;
93 double psat, rhof,rhog;
94 fprops_sat_T(t,&psat,&rhof,&rhog,P,&err);
95 double dt = (Tc-Tt)/NSAT;
96 i = (int)round(((t - Tt)/(Tc - Tt)*(NSAT-1)));
97 double delt = t - ( Tt + i*dt);
98
99 double rhofT,rhogT;
100 rhofT = P->table->satFRho[i] + delt*P->table->satFdRhodt[i];// + 0.5*delt*delt*P->table->satFd2RhodT2[i];
101 rhogT = P->table->satGRho[i] + delt*P->table->satGdRhodt[i];// + 0.5*delt*delt*P->table->satGd2RhodT2[i];
102
103 MSG("SAT TEST ::: %f %f",100*(rhof-rhofT)/rhof,100*(rhog-rhogT)/rhog);
104 */
105 /*****************************************Single Phase Table Testing*****************************************/
106 #define NPOINTS 100000
107 double temp_s = 650;
108 double temp_f = 1650;
109 int nT = NPOINTS;
110 double dT = (temp_f-temp_s)/nT;
111
112
113 double pressH[NPOINTS],enthalpyH[NPOINTS];
114 double entH[NPOINTS], intuH[NPOINTS], gibbsgH[NPOINTS];
115
116 double pressT[NPOINTS],enthalpyT[NPOINTS];
117 double entT[NPOINTS], intuT[NPOINTS], gibbsgT[NPOINTS];
118
119
120 clock_t start = clock();
121 for(i=0; i<nT; ++i){
122 err = FPROPS_NO_ERROR;
123
124 double T = temp_s + i*dT;
125 //FluidState S = fprops_set_Trho(T,rho,P,&err);
126 //pressH[i] = fprops_p(S,&err) ;
127 pressH[i] = P->p_fn(T, rho, P->data,&err) ;
128 enthalpyH[i] = P->h_fn(T, rho, P->data,&err) ;
129 entH[i] = P->s_fn(T, rho, P->data,&err) ;
130 intuH[i] = P->u_fn(T, rho, P->data,&err) ;
131 gibbsgH[i] = P->g_fn(T, rho, P->data,&err) ;
132 }
133
134 clock_t end = clock();
135 double msecH = (double)(end - start) / (CLOCKS_PER_SEC/1000);
136
137 start = clock();
138 for(i=0; i<nT; ++i){
139 err = FPROPS_NO_ERROR;
140
141 double T = temp_s + i*dT;
142
143 pressT[i] = evaluate_ttse_p(P,T, rho) ;
144 enthalpyT[i] = evaluate_ttse_h(P,T, rho) ;
145 entT[i] = evaluate_ttse_s(P,T, rho) ;
146 intuT[i] = evaluate_ttse_u(P,T, rho) ;
147 gibbsgT[i] = evaluate_ttse_g(P,T, rho) ;
148 }
149
150 end = clock();
151 double msecT = (double)(end - start) / (CLOCKS_PER_SEC/1000);
152
153
154 // MSG("Percentage Errors");
155 // MSG("Temp \tPressure \tEnthalpy ");
156 double pererrp,pererrh;
157 for(i=0; i<nT; ++i){
158 //double T = temp_s + i*dT;
159 pererrp = 100*((pressT[i]-pressH[i])/pressH[i]);
160 pererrh = 100*((enthalpyT[i]-enthalpyH[i])/enthalpyH[i]);
161
162 // MSG("%3.6f\t%3.6f\t%3.6f",T, pererrp,pererrh );
163 }
164
165
166 // MSG("Percentage Errors");
167 // MSG("Temp \t\tEntropy \t\tU \t\tG ");
168 double pererrs,pererru,pererrg;
169
170 for(i=0; i<nT; ++i){
171 // double T = temp_s + i*dT;
172 pererrs = 100*((entT[i]-entH[i])/entH[i]);
173 pererru = 100*((intuT[i]-intuH[i])/intuH[i]);
174 pererrg = 100*((gibbsgT[i]-gibbsgH[i])/gibbsgH[i]);
175
176 // MSG("%3.6f\t%3.6f\t%3.6f\t%3.6f",T, pererrs,pererru,pererrg );
177 }
178
179 double av[5]={0,0,0,0,0};
180 i=0;
181 while(i<nT){
182 pererrp = 100*((pressT[i]-pressH[i])/pressH[i]);
183 pererrh = 100*((enthalpyT[i]-enthalpyH[i])/enthalpyH[i]);
184 pererrs = 100*((entT[i]-entH[i])/entH[i]);
185 pererru = 100*((intuT[i]-intuH[i])/intuH[i]);
186 pererrg = 100*((gibbsgT[i]-gibbsgH[i])/gibbsgH[i]);
187 av[0] += fabs(pererrp);
188 av[1] += fabs(pererrh);
189 av[2] += fabs(pererrs);
190 av[3] += fabs(pererru);
191 av[4] += fabs(pererrg);
192 ++i;
193 }
194
195 MSG("AVERAGE percentage errors for 5 variables p h s u and g respectively -->");
196 MSG("%3.6f\t%3.6f\t%3.6f\t%3.6f\t%3.6f",av[0]/nT,av[1]/nT,av[2]/nT,av[3]/nT,av[4]/nT);
197
198 MSG("Helmholtz did %d calculations in %e seconds", nT*5,msecH/1000);
199 MSG("TTSE did %d calculations in %e seconds", nT*5,msecT/1000);
200
201 return 1;
202 }
203
204
205
206
207
208
209
210
211
212

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