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

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

Parent Directory Parent Directory | Revision Log Revision Log


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


1 sid 2982
2     #include <math.h>
3     #include <stdio.h>
4     #include <stdlib.h>
5    
6 sid 2987
7 sid 2982 #include "rundata.h"
8     #include "ttse.h"
9    
10    
11    
12 sid 2987 #define TTSE_DEBUG //sid change
13     #ifdef TTSE_DEBUG
14 sid 2983 # 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 sid 2987 #include <ascend/general/ascMalloc.h>
24     #include <ascend/utilities/error.h>
25 sid 2983
26 sid 2987 #define INIT -23456789.123
27 sid 2982
28    
29 sid 2987 inline TtseMatrix alloc_matrix(int tp, int rhop) {
30 sid 2982
31 sid 2987 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 sid 2982
43    
44 sid 2987 inline void remove_matrix(TtseMatrix mat , int tp){
45 sid 2982
46 sid 2987 ASC_FREE(mat);
47 sid 2982 }
48    
49    
50 sid 2987 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 sid 2982 }
93    
94    
95 sid 2987 void remove_tables(Ttse *table)
96 sid 2982 {
97 sid 2987 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 sid 2989 /*
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 sid 2987
137     void load_tables(PureFluid *P){
138    
139    
140     }
141    
142 sid 2989 /*
143     After building the tables once this should be called to save the files in binary inside tables/
144     */
145 sid 2987 void save_tables(PureFluid *P){
146    
147    
148     }
149    
150 sid 2989 /*
151     Actual building of tables is done here.
152     */
153 sid 2987 void build_tables(PureFluid *P){
154    
155     #ifndef PT
156     #define PT P->table
157    
158 sid 2982 int i,j;
159 sid 2987 FpropsError err = FPROPS_NO_ERROR;
160 sid 2982
161 sid 2987 double tmin,tmax,rhomin,rhomax;
162 sid 2982
163 sid 2989 //Pseudo values for water
164 sid 2990 //Should be implemented else where per fluid
165 sid 2989 PT->tmin = 200;
166     PT->tmax = 800;
167     PT->rhomin = 400;
168     PT->rhomax = 1400;
169    
170 sid 2990
171    
172 sid 2987 tmin = PT->tmin;
173     tmax = PT->tmax;
174     rhomin = PT->rhomin;
175     rhomax = PT->rhomax;
176 sid 2982
177 sid 2987 double dt = (tmax-tmin)/NTP;
178     double drho = (rhomax-rhomin)/NRHOP;
179 sid 2982
180    
181 sid 2992 MSG("DTemp is %f and DRho is %f",dt,drho);
182 sid 2991 MSG("BUILDING TABLES");
183 sid 2982
184 sid 2987 clock_t start = clock();
185 sid 2982
186 sid 2987 for( i = 0; i < NTP; i++)
187     for( j = 0; j < NRHOP; j++){
188 sid 2982
189 sid 2987 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 sid 2989 MSG("Tables built in %f seconds", msec/1000);
220 sid 2987
221    
222    
223     #undef PT
224     #endif
225 sid 2982 }
226    
227    
228 sid 2987
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 sid 2982 {
256    
257 sid 2987 #ifndef PT
258     #define PT P->table
259 sid 2982
260 sid 2987 int i,j;
261     FpropsError err = FPROPS_NO_ERROR;
262 sid 2982
263 sid 2987 double tmin = PT->tmin;
264     double tmax = PT->tmax;
265     double rhomin = PT->rhomin;
266     double rhomax= PT->rhomax;
267 sid 2982
268 sid 2987 double dt = (tmax-tmin)/NTP;
269     double drho = (rhomax-rhomin)/NRHOP;
270 sid 2982
271 sid 2987 i = (int)round(((t-tmin)/(tmax-tmin)*(NTP-1)));
272     j = (int)round(((rho-rhomin)/(rhomax-rhomin)*(NRHOP-1)));
273 sid 2982
274 sid 2987 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 sid 2982 }
293    
294 sid 2987 */
295    
296    
297 sid 2982 void ttse_prepare(PureFluid *P){
298    
299 sid 2987 #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 sid 2983
304    
305 sid 2987 if(!P->table->usettse)
306 sid 2982 return;
307    
308    
309 sid 2987 MSG("Inside TTSE");
310 sid 2982
311 sid 2987 alloc_tables(P->table);
312 sid 2982
313 sid 2987 build_tables(P);
314 sid 2982
315 sid 2987 #ifdef TTSE_DEBUG
316     fclose(F1);
317     #endif
318    
319 sid 2989
320 sid 2982 }
321 sid 2987
322    
323     void ttse_clean(PureFluid *P){
324     remove_tables(P->table);
325 sid 2989
326    
327 sid 2987 }
328    

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