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

Contents of /branches/sid/models/johnpye/fprops/ttse.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2992 - (show annotations) (download) (as text)
Tue Jun 30 06:24:29 2015 UTC (7 years, 5 months ago) by sid
File MIME type: text/x-csrc
File size: 8257 byte(s)


1
2 #include <math.h>
3 #include <stdio.h>
4 #include <stdlib.h>
5
6
7 #include "rundata.h"
8 #include "ttse.h"
9
10
11
12 #define TTSE_DEBUG //sid change
13 #ifdef TTSE_DEBUG
14 # include "color.h"
15 # define MSG FPROPS_MSG
16 # define ERRMSG FPROPS_ERRMSG
17 #else
18 # define MSG(ARGS...) ((void)0)
19 # define ERRMSG(ARGS...) ((void)0)
20 #endif
21
22
23 #include <ascend/general/ascMalloc.h>
24 #include <ascend/utilities/error.h>
25
26 #define INIT -23456789.123
27
28
29 inline TtseMatrix alloc_matrix(int tp, int rhop) {
30
31 TtseMatrix matrix = ASC_NEW_ARRAY (TtseMatrix, tp*rhop); //tp rows and rhop columns
32 /*
33 matrix[10][24] = 99;
34 printf("%f\n", *((double*)matrix+(tp*10)+24));
35 matrix[20][44] = -99;
36 printf("%f\n", *((double*)matrix+(tp*20)+44));
37 matrix[30][64] = 199;
38 printf("%f\n", *((double*)matrix+(tp*30)+64));
39 */
40 return matrix;
41 }
42
43
44 inline void remove_matrix(TtseMatrix mat , int tp){
45
46 ASC_FREE(mat);
47 }
48
49
50 void alloc_tables(Ttse * table)
51 {
52
53 table->s = alloc_matrix(NTP,NRHOP);
54 table->dsdt = alloc_matrix(NTP,NRHOP);
55 table->d2sdt2 = alloc_matrix(NTP,NRHOP);
56 table->dsdrho = alloc_matrix(NTP,NRHOP);
57 table->d2sdrho2 = alloc_matrix(NTP,NRHOP);
58 table->d2sdtdrho = alloc_matrix(NTP,NRHOP);
59
60
61 table->p = alloc_matrix(NTP,NRHOP);
62 table->dpdt = alloc_matrix(NTP,NRHOP);
63 table->d2pdt2 = alloc_matrix(NTP,NRHOP);
64 table->dpdrho = alloc_matrix(NTP,NRHOP);
65 table->d2pdrho2 = alloc_matrix(NTP,NRHOP);
66 table->d2pdtdrho = alloc_matrix(NTP,NRHOP);
67
68
69 table->u = alloc_matrix(NTP,NRHOP);
70 table->dudt = alloc_matrix(NTP,NRHOP);
71 table->d2udt2 = alloc_matrix(NTP,NRHOP);
72 table->dudrho = alloc_matrix(NTP,NRHOP);
73 table->d2udrho2 = alloc_matrix(NTP,NRHOP);
74 table->d2udtdrho = alloc_matrix(NTP,NRHOP);
75
76
77 table->g = alloc_matrix(NTP,NRHOP);
78 table->dgdt = alloc_matrix(NTP,NRHOP);
79 table->d2gdt2 = alloc_matrix(NTP,NRHOP);
80 table->dgdrho = alloc_matrix(NTP,NRHOP);
81 table->d2gdrho2 = alloc_matrix(NTP,NRHOP);
82 table->d2gdtdrho = alloc_matrix(NTP,NRHOP);
83
84
85
86 table->h = alloc_matrix(NTP,NRHOP);
87 table->dhdt = alloc_matrix(NTP,NRHOP);
88 table->d2hdt2 = alloc_matrix(NTP,NRHOP);
89 table->dhdrho = alloc_matrix(NTP,NRHOP);
90 table->d2hdrho2 = alloc_matrix(NTP,NRHOP);
91 table->d2hdtdrho = alloc_matrix(NTP,NRHOP);
92 }
93
94
95 void remove_tables(Ttse *table)
96 {
97 remove_matrix(table->dsdt,NTP);
98 remove_matrix(table->d2sdt2,NTP);
99 remove_matrix(table->dsdrho,NTP);
100 remove_matrix(table->d2sdrho2,NTP);
101 remove_matrix(table->d2sdtdrho,NTP);
102
103 remove_matrix(table->dpdt,NTP);
104 remove_matrix(table->d2pdt2,NTP);
105 remove_matrix(table->dpdrho,NTP);
106 remove_matrix(table->d2pdrho2,NTP);
107 remove_matrix(table->d2pdtdrho,NTP);
108
109
110 remove_matrix(table->dudt,NTP);
111 remove_matrix(table->d2udt2,NTP);
112 remove_matrix(table->dudrho,NTP);
113 remove_matrix(table->d2udrho2,NTP);
114 remove_matrix(table->d2udtdrho,NTP);
115
116 remove_matrix(table->dgdt,NTP);
117 remove_matrix(table->d2gdt2,NTP);
118 remove_matrix(table->dgdrho,NTP);
119 remove_matrix(table->d2gdrho2,NTP);
120 remove_matrix(table->d2gdtdrho,NTP);
121
122 remove_matrix(table->dhdt,NTP);
123 remove_matrix(table->d2hdt2,NTP);
124 remove_matrix(table->dhdrho,NTP);
125 remove_matrix(table->d2hdrho2,NTP);
126 remove_matrix(table->d2hdtdrho,NTP);
127
128 }
129
130
131
132 /*
133 This will load the binary file from tables/ for the liquid of interest and the EOS and populate the matrices.
134 If the files are not present in tables/ then build_tables() should be used.
135 */
136
137 void load_tables(PureFluid *P){
138
139
140 }
141
142 /*
143 After building the tables once this should be called to save the files in binary inside tables/
144 */
145 void save_tables(PureFluid *P){
146
147
148 }
149
150 /*
151 Actual building of tables is done here.
152 */
153 void build_tables(PureFluid *P){
154
155 #ifndef PT
156 #define PT P->table
157
158 int i,j;
159 FpropsError err = FPROPS_NO_ERROR;
160
161 double tmin,tmax,rhomin,rhomax;
162
163 //Pseudo values for water
164 //Should be implemented else where per fluid
165 PT->tmin = 200;
166 PT->tmax = 800;
167 PT->rhomin = 400;
168 PT->rhomax = 1400;
169
170
171
172 tmin = PT->tmin;
173 tmax = PT->tmax;
174 rhomin = PT->rhomin;
175 rhomax = PT->rhomax;
176
177 double dt = (tmax-tmin)/NTP;
178 double drho = (rhomax-rhomin)/NRHOP;
179
180
181 MSG("DTemp is %f and DRho is %f",dt,drho);
182 MSG("BUILDING TABLES");
183
184 clock_t start = clock();
185
186 for( i = 0; i < NTP; i++)
187 for( j = 0; j < NRHOP; j++){
188
189 double t = tmin+i*dt;
190 double rho = rhomin+j*drho;
191
192 PT->p[i][j] = P->p_fn( t, rho , P->data, &err);
193 PT->dpdt[i][j] = P->dpdT_rho_fn( t, rho , P->data, &err);
194 PT->dpdrho[i][j] = P->dpdrho_T_fn( t, rho , P->data, &err);
195 PT->d2pdt2[i][j] = P->d2pdT2_rho_fn( t, rho , P->data, &err);
196 PT->d2pdrho2[i][j] = P->d2pdrho2_T_fn( t, rho , P->data, &err);
197 PT->d2pdtdrho[i][j] = P->d2pdTdrho_fn( t, rho , P->data, &err);
198
199
200 PT->h[i][j] = P->h_fn( t, rho , P->data, &err);
201 PT->dhdt[i][j] = P->dhdT_rho_fn( t, rho , P->data, &err);
202 PT->dhdrho[i][j] = P->dhdrho_T_fn( t, rho , P->data, &err);
203 PT->d2hdt2[i][j] = P->d2hdT2_rho_fn( t, rho , P->data, &err);
204 PT->d2hdrho2[i][j] = P->d2hdrho2_T_fn( t, rho , P->data, &err);
205 PT->d2hdtdrho[i][j] = P->d2hdTdrho_fn( t, rho , P->data, &err);
206
207
208 PT->s[i][j] = P->s_fn( t, rho , P->data, &err);
209 PT->dsdt[i][j] = P->dsdT_rho_fn( t, rho , P->data, &err);
210 PT->dsdrho[i][j] = P->dsdrho_T_fn( t, rho , P->data, &err);
211 PT->d2sdt2[i][j] = P->d2sdT2_rho_fn( t, rho , P->data, &err);
212 PT->d2sdrho2[i][j] = P->d2sdrho2_T_fn( t, rho , P->data, &err);
213 PT->d2sdtdrho[i][j] = P->d2sdTdrho_fn( t, rho , P->data, &err);
214 }
215
216
217 clock_t end = clock();
218 double msec = (double)(end - start) / (CLOCKS_PER_SEC/1000);
219 MSG("Tables built in %f seconds", msec/1000);
220
221
222
223 #undef PT
224 #endif
225 }
226
227
228
229 #define EVALTTSEFN(VAR) \
230 double evaluate_ttse_##VAR(PureFluid *P , double t, double rho){\
231 int i,j;\
232 double tmin = P->table->tmin; double tmax = P->table->tmax;\
233 double rhomin = P->table->rhomin; double rhomax= P->table->rhomax;\
234 double dt = (tmax-tmin)/NTP;\
235 double drho = (rhomax-rhomin)/NRHOP;\
236 i = (int)round(((t-tmin)/(tmax-tmin)*(NTP-1)));\
237 j = (int)round(((rho-rhomin)/(rhomax-rhomin)*(NRHOP-1)));\
238 double delt = t - ( tmin + i*dt);\
239 double delrho = rho - ( rhomin + j*drho);\
240 double ttse##VAR = P->table->VAR[i][j]+ 0.5*delt*delt*P->table->d2##VAR##dt2[i][j]\
241 + delrho*P->table->d##VAR##drho[i][j] + 0.5*delrho*delrho*P->table->d2##VAR##drho2[i][j]\
242 + delrho*delt*P->table->d2##VAR##dtdrho[i][j];\
243 return ttse##VAR;\
244 }
245
246
247 EVALTTSEFN(p);
248 EVALTTSEFN(h);
249 EVALTTSEFN(s);
250
251
252
253 /*
254 double evaluate_ttse_p(PureFluid *P , double t, double rho)
255 {
256
257 #ifndef PT
258 #define PT P->table
259
260 int i,j;
261 FpropsError err = FPROPS_NO_ERROR;
262
263 double tmin = PT->tmin;
264 double tmax = PT->tmax;
265 double rhomin = PT->rhomin;
266 double rhomax= PT->rhomax;
267
268 double dt = (tmax-tmin)/NTP;
269 double drho = (rhomax-rhomin)/NRHOP;
270
271 i = (int)round(((t-tmin)/(tmax-tmin)*(NTP-1)));
272 j = (int)round(((rho-rhomin)/(rhomax-rhomin)*(NRHOP-1)));
273
274 double delt = t - ( tmin + i*dt);
275 double delrho = rho - ( rhomin + j*drho);
276
277 MSG("%f %f",i ,j ,delt ,delrho );
278 MSG("%f %f %f %f %f %f", PT->p[i][j], PT->dpdt[i][j], PT->d2pdt2[i][j], PT->dpdrho[i][j], PT->d2pdrho2[i][j], PT->d2pdtdrho[i][j]);
279
280
281 double ttse_p = PT->p[i][j] + delt*PT->dpdt[i][j] + 0.5*delt*delt*PT->d2pdt2[i][j]
282 + delrho*PT->dpdrho[i][j] + 0.5*delrho*delrho*PT->d2pdrho2[i][j]
283 + delrho*delt*PT->d2pdtdrho[i][j];
284
285
286 MSG("%e %e", ttse_p , P->p_fn(t, rho, P->data,&err) );
287
288
289 return(ttse_p)
290 #undef PT
291 #endif
292 }
293
294 */
295
296
297 void ttse_prepare(PureFluid *P){
298
299 #ifdef TTSE_DEBUG
300 FILE *F1 = fopen("ttse.txt","w");
301 //fprintf(F1,"%f %f\n",t, P->p_fn( t, rho , P->data,&err) );
302 #endif
303
304
305 if(!P->table->usettse)
306 return;
307
308
309 MSG("Inside TTSE");
310
311 alloc_tables(P->table);
312
313 build_tables(P);
314
315 #ifdef TTSE_DEBUG
316 fclose(F1);
317 #endif
318
319
320 }
321
322
323 void ttse_clean(PureFluid *P){
324 remove_tables(P->table);
325
326
327 }
328

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