johnpye/fprops/ttse.c


#include <math.h>
#include <stdio.h>
#include <stdlib.h>


#include "rundata.h"
#include "ttse.h"


#define TTSE_DEBUG //sid change
#ifdef TTSE_DEBUG
# include "color.h"
# define MSG FPROPS_MSG
# define ERRMSG FPROPS_ERRMSG
#else
# define MSG(ARGS...) ((void)0)
# define ERRMSG(ARGS...) ((void)0)
#endif


#include <ascend/general/ascMalloc.h>
#include <ascend/utilities/error.h>

#define INIT -23456789.123


inline TtseMatrix alloc_matrix(int tp, int rhop) {

    TtseMatrix matrix =  ASC_NEW_ARRAY (TtseMatrix, tp*rhop);  //tp rows and rhop columns
/*
    matrix[10][24] = 99;
    printf("%f\n", *((double*)matrix+(tp*10)+24));
    matrix[20][44] = -99;
    printf("%f\n", *((double*)matrix+(tp*20)+44));
    matrix[30][64] = 199;
    printf("%f\n", *((double*)matrix+(tp*30)+64));
*/
    return matrix;
}


inline void remove_matrix(TtseMatrix mat , int tp){

    ASC_FREE(mat);
}


void alloc_tables(Ttse * table)
{

    table->s = alloc_matrix(NTP,NRHOP);
    table->dsdt = alloc_matrix(NTP,NRHOP);
    table->d2sdt2 = alloc_matrix(NTP,NRHOP);
    table->dsdrho = alloc_matrix(NTP,NRHOP);
    table->d2sdrho2 = alloc_matrix(NTP,NRHOP);
    table->d2sdtdrho = alloc_matrix(NTP,NRHOP);


    table->p = alloc_matrix(NTP,NRHOP);
    table->dpdt = alloc_matrix(NTP,NRHOP);
    table->d2pdt2 = alloc_matrix(NTP,NRHOP);
    table->dpdrho = alloc_matrix(NTP,NRHOP);
    table->d2pdrho2 = alloc_matrix(NTP,NRHOP);
    table->d2pdtdrho = alloc_matrix(NTP,NRHOP);


    table->u = alloc_matrix(NTP,NRHOP);
    table->dudt = alloc_matrix(NTP,NRHOP);
    table->d2udt2 = alloc_matrix(NTP,NRHOP);
    table->dudrho = alloc_matrix(NTP,NRHOP);
    table->d2udrho2 = alloc_matrix(NTP,NRHOP);
    table->d2udtdrho = alloc_matrix(NTP,NRHOP);


    table->g = alloc_matrix(NTP,NRHOP);
    table->dgdt = alloc_matrix(NTP,NRHOP);
    table->d2gdt2 = alloc_matrix(NTP,NRHOP);
    table->dgdrho = alloc_matrix(NTP,NRHOP);
    table->d2gdrho2 = alloc_matrix(NTP,NRHOP);
    table->d2gdtdrho = alloc_matrix(NTP,NRHOP);


    table->h = alloc_matrix(NTP,NRHOP);
    table->dhdt = alloc_matrix(NTP,NRHOP);
    table->d2hdt2 = alloc_matrix(NTP,NRHOP);
    table->dhdrho = alloc_matrix(NTP,NRHOP);
    table->d2hdrho2 = alloc_matrix(NTP,NRHOP);
    table->d2hdtdrho = alloc_matrix(NTP,NRHOP);
}


void remove_tables(Ttse *table)
{
    remove_matrix(table->dsdt,NTP);
    remove_matrix(table->d2sdt2,NTP);
    remove_matrix(table->dsdrho,NTP);
    remove_matrix(table->d2sdrho2,NTP);
    remove_matrix(table->d2sdtdrho,NTP);

    remove_matrix(table->dpdt,NTP);
    remove_matrix(table->d2pdt2,NTP);
    remove_matrix(table->dpdrho,NTP);
    remove_matrix(table->d2pdrho2,NTP);
    remove_matrix(table->d2pdtdrho,NTP);


    remove_matrix(table->dudt,NTP);
    remove_matrix(table->d2udt2,NTP);
    remove_matrix(table->dudrho,NTP);
    remove_matrix(table->d2udrho2,NTP);
    remove_matrix(table->d2udtdrho,NTP);

    remove_matrix(table->dgdt,NTP);
    remove_matrix(table->d2gdt2,NTP);
    remove_matrix(table->dgdrho,NTP);
    remove_matrix(table->d2gdrho2,NTP);
    remove_matrix(table->d2gdtdrho,NTP);

    remove_matrix(table->dhdt,NTP);
    remove_matrix(table->d2hdt2,NTP);
    remove_matrix(table->dhdrho,NTP);
    remove_matrix(table->d2hdrho2,NTP);
    remove_matrix(table->d2hdtdrho,NTP);

}


/*
    This will load the binary file from tables/ for the liquid of interest and the EOS and populate the matrices.
    If the files are not present in tables/ then build_tables() should be used.
*/

void load_tables(PureFluid *P){


}

/*
    After building the tables once this should be called to save the files in binary inside tables/
*/
void save_tables(PureFluid *P){


}

/*
    Actual building of tables is done here.
*/
void build_tables(PureFluid *P){

    #ifndef PT
    #define PT P->table

    int i,j;
    FpropsError err = FPROPS_NO_ERROR;

    double tmin,tmax,rhomin,rhomax;

//Pseudo values for water
//Should be implemented else where per fluid
    PT->tmin = 200;
    PT->tmax = 800;
    PT->rhomin = 400;
    PT->rhomax = 1400;


    tmin = PT->tmin;
    tmax = PT->tmax;
    rhomin = PT->rhomin;
    rhomax = PT->rhomax;

    double dt = (tmax-tmin)/NTP;
    double drho = (rhomax-rhomin)/NRHOP;


    MSG("DTemp is %f and DRho is  %f",dt,drho);
    MSG("BUILDING TABLES");

    clock_t start = clock();

    for( i = 0; i < NTP; i++)
    for( j = 0; j < NRHOP; j++){

        double t  = tmin+i*dt;
        double rho  = rhomin+j*drho;

        PT->p[i][j] = P->p_fn( t, rho , P->data, &err);
        PT->dpdt[i][j] = P->dpdT_rho_fn( t, rho , P->data, &err);
        PT->dpdrho[i][j] = P->dpdrho_T_fn( t, rho , P->data, &err);
        PT->d2pdt2[i][j] = P->d2pdT2_rho_fn( t, rho , P->data, &err);
        PT->d2pdrho2[i][j] = P->d2pdrho2_T_fn( t, rho , P->data, &err);
        PT->d2pdtdrho[i][j] = P->d2pdTdrho_fn( t, rho , P->data, &err);


        PT->h[i][j] = P->h_fn( t, rho , P->data, &err);
        PT->dhdt[i][j] = P->dhdT_rho_fn( t, rho , P->data, &err);
        PT->dhdrho[i][j] = P->dhdrho_T_fn( t, rho , P->data, &err);
        PT->d2hdt2[i][j] = P->d2hdT2_rho_fn( t, rho , P->data, &err);
        PT->d2hdrho2[i][j] = P->d2hdrho2_T_fn( t, rho , P->data, &err);
        PT->d2hdtdrho[i][j] = P->d2hdTdrho_fn( t, rho , P->data, &err);


        PT->s[i][j] = P->s_fn( t, rho , P->data, &err);
        PT->dsdt[i][j] = P->dsdT_rho_fn( t, rho , P->data, &err);
        PT->dsdrho[i][j] = P->dsdrho_T_fn( t, rho , P->data, &err);
        PT->d2sdt2[i][j] = P->d2sdT2_rho_fn( t, rho , P->data, &err);
        PT->d2sdrho2[i][j] = P->d2sdrho2_T_fn( t, rho , P->data, &err);
        PT->d2sdtdrho[i][j] = P->d2sdTdrho_fn( t, rho , P->data, &err);
    }


    clock_t end = clock();
    double msec = (double)(end - start) / (CLOCKS_PER_SEC/1000);
    MSG("Tables built in %f seconds", msec/1000);


    #undef PT
    #endif
}


#define EVALTTSEFN(VAR) \
    double evaluate_ttse_##VAR(PureFluid *P , double t, double rho){\
            int i,j;\
            double tmin = P->table->tmin; double tmax = P->table->tmax;\
            double rhomin  = P->table->rhomin; double rhomax= P->table->rhomax;\
            double dt = (tmax-tmin)/NTP;\
            double drho = (rhomax-rhomin)/NRHOP;\
            i = (int)round(((t-tmin)/(tmax-tmin)*(NTP-1)));\
            j = (int)round(((rho-rhomin)/(rhomax-rhomin)*(NRHOP-1)));\
            double delt = t - ( tmin + i*dt);\
            double delrho = rho - ( rhomin + j*drho);\
            double ttse##VAR = P->table->VAR[i][j]+ 0.5*delt*delt*P->table->d2##VAR##dt2[i][j]\
                 + delrho*P->table->d##VAR##drho[i][j] + 0.5*delrho*delrho*P->table->d2##VAR##drho2[i][j]\
                 + delrho*delt*P->table->d2##VAR##dtdrho[i][j];\
            return ttse##VAR;\
        }


EVALTTSEFN(p);
EVALTTSEFN(h);
EVALTTSEFN(s);


/*
double evaluate_ttse_p(PureFluid *P , double t, double rho)
{

    #ifndef PT
    #define PT P->table

    int i,j;
    FpropsError err = FPROPS_NO_ERROR;

    double tmin = PT->tmin;
    double tmax = PT->tmax;
    double rhomin  = PT->rhomin;
    double rhomax= PT->rhomax;

    double dt = (tmax-tmin)/NTP;
    double drho = (rhomax-rhomin)/NRHOP;

    i = (int)round(((t-tmin)/(tmax-tmin)*(NTP-1)));
    j = (int)round(((rho-rhomin)/(rhomax-rhomin)*(NRHOP-1)));

    double delt = t - ( tmin + i*dt);
    double delrho = rho - ( rhomin + j*drho);

    MSG("%f %f",i ,j ,delt ,delrho );
    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]);


    double ttse_p = PT->p[i][j] + delt*PT->dpdt[i][j] + 0.5*delt*delt*PT->d2pdt2[i][j]
            + delrho*PT->dpdrho[i][j] + 0.5*delrho*delrho*PT->d2pdrho2[i][j]
            + delrho*delt*PT->d2pdtdrho[i][j];


    MSG("%e  %e", ttse_p , P->p_fn(t, rho, P->data,&err)  );


    return(ttse_p)
    #undef PT
    #endif
}

*/


void ttse_prepare(PureFluid *P){

#ifdef TTSE_DEBUG
    FILE *F1 = fopen("ttse.txt","w");
    //fprintf(F1,"%f   %f\n",t, P->p_fn( t, rho , P->data,&err) );
#endif


    if(!P->table->usettse)
        return;


    MSG("Inside TTSE");

    alloc_tables(P->table);

    build_tables(P);

#ifdef TTSE_DEBUG
    fclose(F1);
#endif


}


void ttse_clean(PureFluid *P){
    remove_tables(P->table);


}

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.5deltdelt*P->table->d2##VAR##dt2[i][j]\
241	+ delrhoP->table->d##VAR##drho[i][j] + 0.5delrhodelrhoP->table->d2##VAR##drho2[i][j]\
242	+ delrhodeltP->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] + deltPT->dpdt[i][j] + 0.5deltdeltPT->d2pdt2[i][j]
282	+ delrhoPT->dpdrho[i][j] + 0.5delrhodelrhoPT->d2pdrho2[i][j]
283	+ delrhodeltPT->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