/[ascend]/branches/vishnu/models/johnpye/fprops/vishnu/scratch/main.c
ViewVC logotype

Annotation of /branches/vishnu/models/johnpye/fprops/vishnu/scratch/main.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3168 - (hide annotations) (download) (as text)
Wed Aug 17 06:29:34 2016 UTC (5 years, 11 months ago) by vishnu
File MIME type: text/x-csrc
File size: 7820 byte(s)
pre-final
1 vishnu 3115
2     #include <stdio.h>
3     #include <stdlib.h>
4     #include <string.h>
5     #include <assert.h>
6 vishnu 3116 #include <math.h>
7 vishnu 3115
8     #define BUF 1024
9    
10     #define MAXC 10
11    
12 vishnu 3120 #define PATH "incomp_liq_data/LiBr.dat"
13 vishnu 3115
14     #define ABORT \
15 vishnu 3118 printf("\nExiting Program!\n"); \
16 vishnu 3115 exit(0);
17    
18 vishnu 3120 // json parser
19 vishnu 3115 #define READ(Q) \
20     pos = ftell(in); \
21     while(!strstr(word,"},")) { \
22     fgets(word,BUF,in); \
23     if(strstr(word,"\"type\": \"notdefined\"")) { \
24     test_liq->Q.type = (char*)malloc(sizeof("notdefined")); \
25     strcpy(test_liq->Q.type,"notdefined"); \
26     } \
27     else if(strstr(word,"\"type\": \"polynomial\"")) { \
28     test_liq->Q.type = (char*)malloc(sizeof("polynomial")); \
29     strcpy(test_liq->Q.type,"polynomial"); \
30     } \
31     else if(strstr(word,"\"type\": \"exppolynomial\"")) { \
32     test_liq->Q.type = (char*)malloc(sizeof("exppolynomial")); \
33     strcpy(test_liq->Q.type,"exppolynomial"); \
34     } \
35 vishnu 3116 else if(strstr(word,"\"type\": \"exponential\"")) { \
36     test_liq->Q.type = (char*)malloc(sizeof("exponential")); \
37     strcpy(test_liq->Q.type,"exponential"); \
38     } \
39 vishnu 3115 } \
40 vishnu 3120 if(strcmp(test_liq->Q.type,"notdefined")) { \
41 vishnu 3115 fseek(in,pos,SEEK_SET); \
42     fgets(word,BUF,in); \
43 vishnu 3120 double coefs[MAXC][MAXC]; \
44     int numc_row = 0, numc_col = 0; \
45 vishnu 3115 while(!strstr(word,"},")) { \
46 vishnu 3120 if((strstr(word,"e+")||strstr(word,"e-"))&&!strstr(word,"NRMS")) { \
47     coefs[numc_row][numc_col] = to_number(word); \
48     numc_col++; \
49     assert(numc_col<MAXC); \
50 vishnu 3115 } \
51 vishnu 3120 else if(strstr(word,"],")) { \
52     numc_row++; \
53     test_liq->Q.numc_c = numc_col; \
54     numc_col = 0; \
55     assert(numc_row<MAXC); \
56 vishnu 3115 } \
57     fgets(word,BUF,in); \
58     } \
59 vishnu 3120 test_liq->Q.numc_r = numc_row; \
60     test_liq->Q.coeff = (double**)malloc(test_liq->Q.numc_r*sizeof(double*)); \
61     int i,j; \
62     for(i=0;i<test_liq->Q.numc_r;i++) \
63     test_liq->Q.coeff[i] = (double*)malloc(test_liq->Q.numc_c*sizeof(double)); \
64     for(i=0;i<test_liq->Q.numc_r;i++) \
65     for(j=0;j<test_liq->Q.numc_c;j++) \
66     test_liq->Q.coeff[i][j] = coefs[i][j]; \
67 vishnu 3115 }
68    
69     typedef struct {
70    
71     char* type;
72 vishnu 3120 double** coeff;
73     int numc_r, numc_c; // number of rows and columns in coefficient matrix respectively
74 vishnu 3115
75     } coefficients;
76    
77     typedef struct {
78    
79 vishnu 3120 coefficients T_freeze, conductivity, density, specific_heat, viscosity, saturation_pressure;
80 vishnu 3115 char* x_id;
81     char* description;
82     double T_base, T_max, T_min, T_minPsat, x_base, x_max, x_min;
83    
84     } fprops;
85    
86     double extracted_number(char w[BUF]) {
87    
88     char *token;
89    
90     /* get the first token */
91     token = strtok(w, " ,");
92    
93     /* get the second token which is the number*/
94     token = strtok(NULL, " ,");
95    
96     return (atof(token));
97    
98     }
99    
100     double to_number(char w[BUF]) {
101    
102     char *token;
103    
104     /* get the first token */
105 vishnu 3118 token = strtok(w, " ,");
106 vishnu 3115
107     return (atof(token));
108    
109 vishnu 3116 }
110    
111 vishnu 3120 // evaluation function for polynamial type coefficients
112 vishnu 3116 double eval_poly(coefficients c, double T, double x) {
113    
114     double res = 0.0;
115 vishnu 3120 int i,j;
116 vishnu 3116
117 vishnu 3120 for(i=0;i<c.numc_r;i++)
118     for(j=0;j<c.numc_c;j++)
119     res += pow(x,j)*c.coeff[i][j]*pow(T,i);
120 vishnu 3116
121     return res;
122    
123     }
124    
125 vishnu 3120 // evaluation function for exponential-polynomial type coefficients
126 vishnu 3116 double eval_exppoly(coefficients c, double T, double x) {
127    
128     return exp(eval_poly(c,T,x));
129    
130     }
131    
132 vishnu 3120 // evaluation function for exponential type coefficients
133 vishnu 3116 double eval_expo(coefficients c, double T, double x) {
134    
135 vishnu 3120 return exp(c.coeff[0][0]/(T+c.coeff[0][1])-c.coeff[0][2]);
136 vishnu 3116
137 vishnu 3115 }
138    
139 vishnu 3116 // generalized property evaluation functions: macro definition
140     #define PROP_EVAL(Q) \
141     double eval_ ## Q (fprops *fluid, double T, double x) { \
142     assert(T>0&&T>=fluid->T_min&&T<=fluid->T_max); \
143     assert(x>=fluid->x_min&&x<=fluid->x_max); \
144 vishnu 3118 if(!strcmp(fluid->Q.type,"polynomial")) return eval_poly(fluid->Q,T-fluid->T_base,x-fluid->x_base); \
145     else if(!strcmp(fluid->Q.type,"exppolynomial")) return eval_exppoly(fluid->Q,T-fluid->T_base,x-fluid->x_base); \
146 vishnu 3120 else if(!strcmp(fluid->Q.type,"exponential")) return eval_expo(fluid->Q,T,x); \
147 vishnu 3116 else { \
148     printf("\nType not defined.\n"); \
149     ABORT \
150     } \
151     }
152    
153     // declaration and definition of property evaluation functions for conductivity, density, specific-heat and viscosity
154     PROP_EVAL(conductivity)
155     PROP_EVAL(density)
156     PROP_EVAL(specific_heat)
157     PROP_EVAL(viscosity)
158 vishnu 3120 PROP_EVAL(saturation_pressure)
159 vishnu 3116
160 vishnu 3168 // All incompressible fluids have an arbitrary reference state for enthalpy and entropy. During initialisation, the reference state is defined as a temperature of 20 ��C and a pressure of 1 atm according to the U.S. National Institute of Standards and Technology.
161    
162     double eval_u(fprops *fluid, double T, double x) {
163     assert(!strcmp(fluid->specific_heat.type,"polynomial"));
164     double Tref = 293.15; // 20 ��C
165     // We dont know which if the three statements below are required... no validation data as of now
166     //T = T-fluid->T_base;
167     //Tref = Tref-fluid->T_base;
168     //x = x-fluid->x_base;
169     double u = 0.0;
170     int i,j;
171     for(i=0;i<fluid->specific_heat.numc_r;i++)
172     for(j=0;j<fluid->specific_heat.numc_c;j++)
173     u += pow(x,j)/(i+1.0)*fluid->specific_heat.coeff[i][j]*(pow(T,i+1)-pow(Tref,i+1));
174     return u;
175    
176     }
177    
178     double eval_s(fprops *fluid, double T, double x) {
179     assert(!strcmp(fluid->specific_heat.type,"polynomial"));
180     double Tref = 293.15; // 20 ��C
181     // We dont know which if the three statements below are required... no validation data as of now
182     //T = T-fluid->T_base;
183     //Tref = Tref-fluid->T_base;
184     //x = x-fluid->x_base;
185     double s = 0.0;
186     int i,j;
187     double log_term = log(T/Tref);
188     for(i=0;i<fluid->specific_heat.numc_r-1;i++)
189     for(j=0;j<fluid->specific_heat.numc_c;j++)
190     s += pow(x,j)*(fluid->specific_heat.coeff[0][j]*log_term+1.0/(i+1.0)*fluid->specific_heat.coeff[i+1][j]*(pow(T,i+1)-pow(Tref,i+1)));
191     return s;
192    
193     }
194    
195 vishnu 3115 int main() {
196    
197 vishnu 3116 // parsing from input (json format) and loading data to new fluid data structures
198 vishnu 3115 FILE* in;
199    
200     fprops* test_liq;
201     test_liq = (fprops*)malloc(sizeof(fprops));
202    
203     long int pos;
204    
205     if(!(in = fopen(PATH,"r"))) {
206    
207     printf("\nInput file cannot be located.\n");
208     ABORT
209    
210     }
211    
212     char word[BUF];
213    
214     while(!feof(in)) {
215    
216 vishnu 3116 fgets(word,BUF,in);
217 vishnu 3115
218 vishnu 3116 if(strstr(word,"\"T_freeze\"")) { READ(T_freeze); } // not required for incompressible EOS; implemented to check input coefficients of "type": "notdefined"
219     else if(strstr(word,"\"Tbase\"")) test_liq->T_base = extracted_number(word);
220     else if(strstr(word,"\"Tmax\"")) test_liq->T_max = extracted_number(word);
221     else if(strstr(word,"\"Tmin\"")) test_liq->T_min = extracted_number(word);
222     else if(strstr(word,"\"TminPsat\"")) test_liq->T_minPsat = extracted_number(word);
223     else if(strstr(word,"\"conductivity\"")) { READ(conductivity); }
224     else if(strstr(word,"\"density\"")) { READ(density); }
225     else if(strstr(word,"\"specific_heat\"")) { READ(specific_heat); }
226     else if(strstr(word,"\"viscosity\"")) { READ(viscosity); }
227 vishnu 3120 else if(strstr(word,"\"saturation_pressure\"")) { READ(saturation_pressure); }
228 vishnu 3116 else if(strstr(word,"\"xbase\"")) test_liq->x_base = extracted_number(word);
229     else if(strstr(word,"\"xmax\"")) test_liq->x_max = extracted_number(word);
230     else if(strstr(word,"\"xmin\"")) test_liq->x_min = extracted_number(word);
231 vishnu 3115 else if(strstr(word,"\"xid\": \"pure\"")) {
232     test_liq->x_id = (char*)malloc(sizeof("pure"));
233     strcpy(test_liq->x_id,"pure");
234     }
235    
236     }
237    
238     fclose(in);
239    
240 vishnu 3116 // test
241     double T = test_liq->T_min;
242    
243     FILE *out;
244     out = fopen("test_res.dat","w");
245    
246 vishnu 3168 fprintf(out,"VARIABLES = \"Temperature [C]\", \"Conductivity [W/m/K]\", \"Density [kg/m<sup>3</sup>]\", \"Heat Capacity [J/Kg/K]\", \"Viscosity [Pa s]\", \"Saturation Pressure [Pa]\", \"Internal Energy [SI]\", \"Entropy [SI]\"\n");
247 vishnu 3116 do {
248    
249 vishnu 3168 fprintf(out,"%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\n",T-273.15,eval_conductivity(test_liq, T, 0.39),eval_density(test_liq, T, 0.39),eval_specific_heat(test_liq, T, 0.39),eval_viscosity(test_liq, T, 0.39),eval_saturation_pressure(test_liq, T, 0.39),eval_u(test_liq, T, 0.39),eval_s(test_liq, T, 0.39));
250 vishnu 3116 T+=1;
251    
252     } while(T<=test_liq->T_max);
253    
254     fclose(out);
255    
256 vishnu 3115 return 0;
257    
258     }

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