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	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.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	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] + 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	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