/[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 3024 - (show annotations) (download) (as text)
Thu Jul 23 04:16:35 2015 UTC (3 years, 11 months ago) by sid
File MIME type: text/x-csrc
File size: 5539 byte(s)
first attempt at two phase

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 #define SPOINTS 1000
59
60 double Tt = P->data->T_t;
61 double Tc = P->data->T_c;
62 double dt2p = (Tc - Tt)/SPOINTS;
63
64 printf("{");
65 for(i=0;i<SPOINTS+1;i++)
66 {
67 double T = Tt + i*dt2p;
68 double psat, rhof,rhog;
69 fprops_sat_T(T,&psat,&rhof,&rhog,P,&err);
70
71 if(i!= SPOINTS)
72 printf("{%f, %f, %f},\n",T, rhof,rhog);
73 else
74 printf("{%f, %f, %f}};\n",T, rhof,rhog);
75 }
76
77 exit(1);
78 */
79
80
81
82
83
84 /*****************************************Two Phase Table Testing*****************************************/
85 double Tt = P->data->T_t;
86 double Tc = P->data->T_c;
87 double t = Tc - 53.341233;
88 double psat, rhof,rhog;
89 fprops_sat_T(t,&psat,&rhof,&rhog,P,&err);
90 double dT = (Tc-Tt)/NSAT;
91 i = (int)round(((t - Tt)/(Tc - Tt)*(NSAT-1)));
92 double delt = t - ( Tt + i*dT);
93
94 double rhofT,rhogT;
95 rhofT = P->table->satFRho[i] + delt*P->table->satFdRhodt[i];// + 0.5*delt*delt*P->table->satFd2RhodT2[i];
96 rhogT = P->table->satGRho[i] + delt*P->table->satGdRhodt[i];// + 0.5*delt*delt*P->table->satGd2RhodT2[i];
97
98 MSG("SAT TEST ::: %f %f",100*(rhof-rhofT)/rhof,100*(rhog-rhogT)/rhog);
99
100 /*****************************************Single Phase Table Testing*****************************************/
101 #define NPOINTS 100000
102 double temp_s = 650;
103 double temp_f = 1650;
104 int nT = NPOINTS;
105 dT = (temp_f-temp_s)/nT;
106
107
108 double pressH[NPOINTS],enthalpyH[NPOINTS];
109 double entH[NPOINTS], intuH[NPOINTS], gibbsgH[NPOINTS];
110
111 double pressT[NPOINTS],enthalpyT[NPOINTS];
112 double entT[NPOINTS], intuT[NPOINTS], gibbsgT[NPOINTS];
113
114
115 clock_t start = clock();
116 for(i=0; i<nT; ++i){
117 err = FPROPS_NO_ERROR;
118
119 double T = temp_s + i*dT;
120 //FluidState S = fprops_set_Trho(T,rho,P,&err);
121 //pressH[i] = fprops_p(S,&err) ;
122 pressH[i] = P->p_fn(T, rho, P->data,&err) ;
123 enthalpyH[i] = P->h_fn(T, rho, P->data,&err) ;
124 entH[i] = P->s_fn(T, rho, P->data,&err) ;
125 intuH[i] = P->u_fn(T, rho, P->data,&err) ;
126 gibbsgH[i] = P->g_fn(T, rho, P->data,&err) ;
127 }
128
129 clock_t end = clock();
130 double msecH = (double)(end - start) / (CLOCKS_PER_SEC/1000);
131
132 start = clock();
133 for(i=0; i<nT; ++i){
134 err = FPROPS_NO_ERROR;
135
136 double T = temp_s + i*dT;
137
138 pressT[i] = evaluate_ttse_p(P,T, rho) ;
139 enthalpyT[i] = evaluate_ttse_h(P,T, rho) ;
140 entT[i] = evaluate_ttse_s(P,T, rho) ;
141 intuT[i] = evaluate_ttse_u(P,T, rho) ;
142 gibbsgT[i] = evaluate_ttse_g(P,T, rho) ;
143 }
144
145 end = clock();
146 double msecT = (double)(end - start) / (CLOCKS_PER_SEC/1000);
147
148
149 // MSG("Percentage Errors");
150 // MSG("Temp \tPressure \tEnthalpy ");
151 double pererrp,pererrh;
152 for(i=0; i<nT; ++i){
153 //double T = temp_s + i*dT;
154 pererrp = 100*((pressT[i]-pressH[i])/pressH[i]);
155 pererrh = 100*((enthalpyT[i]-enthalpyH[i])/enthalpyH[i]);
156
157 // MSG("%3.6f\t%3.6f\t%3.6f",T, pererrp,pererrh );
158 }
159
160
161 // MSG("Percentage Errors");
162 // MSG("Temp \t\tEntropy \t\tU \t\tG ");
163 double pererrs,pererru,pererrg;
164
165 for(i=0; i<nT; ++i){
166 // double T = temp_s + i*dT;
167 pererrs = 100*((entT[i]-entH[i])/entH[i]);
168 pererru = 100*((intuT[i]-intuH[i])/intuH[i]);
169 pererrg = 100*((gibbsgT[i]-gibbsgH[i])/gibbsgH[i]);
170
171 // MSG("%3.6f\t%3.6f\t%3.6f\t%3.6f",T, pererrs,pererru,pererrg );
172 }
173
174 double av[5]={0,0,0,0,0};
175 i=0;
176 while(i<nT){
177 pererrp = 100*((pressT[i]-pressH[i])/pressH[i]);
178 pererrh = 100*((enthalpyT[i]-enthalpyH[i])/enthalpyH[i]);
179 pererrs = 100*((entT[i]-entH[i])/entH[i]);
180 pererru = 100*((intuT[i]-intuH[i])/intuH[i]);
181 pererrg = 100*((gibbsgT[i]-gibbsgH[i])/gibbsgH[i]);
182 av[0] += fabs(pererrp);
183 av[1] += fabs(pererrh);
184 av[2] += fabs(pererrs);
185 av[3] += fabs(pererru);
186 av[4] += fabs(pererrg);
187 ++i;
188 }
189
190 MSG("AVERAGE percentage errors for 5 variables p h s u and g respectively -->");
191 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);
192
193 MSG("Helmholtz did %d calculations in %e seconds", nT*5,msecH/1000);
194 MSG("TTSE did %d calculations in %e seconds", nT*5,msecT/1000);
195
196 return 1;
197 }
198
199
200
201
202
203
204
205
206
207

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