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

Diff of /branches/sid/models/johnpye/fprops/test/test_ttse.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3023 by sid, Thu Jul 2 01:47:21 2015 UTC revision 3024 by sid, Thu Jul 23 04:16:35 2015 UTC
# Line 34  Line 34 
34  #define TOL_RHO 1e-3  #define TOL_RHO 1e-3
35    
36  int main(void){  int main(void){
37       PureFluid *P;      PureFluid *P;
38      FpropsError err;      FpropsError err;
39      const char *helmfluids[] = { "water"};      const char *helmfluids[] = { "water"};
40      //const int n = sizeof(helmfluids)/sizeof(char *);      //const int n = sizeof(helmfluids)/sizeof(char *);
# Line 45  int main(void){ Line 45  int main(void){
45    
46      P = (PureFluid *)fprops_fluid(helmfluids[0],"helmholtz",NULL);      P = (PureFluid *)fprops_fluid(helmfluids[0],"helmholtz",NULL);
47    
48    
49        MSG("Triple Point Tt %f  & Crit Temp  %f",P->data->T_t,P->data->T_c);
50    
51    
52      MSG("Comparing Helmholtz vs TTSE");      MSG("Comparing Helmholtz vs TTSE");
53    
54    
55      double rho = 1200;      double rho = 1000;
56    
57    /* //Plot in Mathematica Saturation dome
58        #define SPOINTS 1000
59    
60        double Tt = P->data->T_t;
61        double Tc = P->data->T_c;
62        double dt2p = (Tc - Tt)/SPOINTS;
63    
64        printf("{");
65        for(i=0;i<SPOINTS+1;i++)
66        {
67            double T = Tt + i*dt2p;
68            double psat, rhof,rhog;
69            fprops_sat_T(T,&psat,&rhof,&rhog,P,&err);
70    
71            if(i!= SPOINTS)
72                printf("{%f, %f, %f},\n",T, rhof,rhog);
73            else
74                printf("{%f, %f, %f}};\n",T, rhof,rhog);
75        }
76    
77        exit(1);
78    */
79    
80    
81    
82    
83    
84    /*****************************************Two Phase Table Testing*****************************************/
85        double Tt = P->data->T_t;
86        double Tc = P->data->T_c;
87        double t = Tc - 53.341233;
88        double psat, rhof,rhog;
89        fprops_sat_T(t,&psat,&rhof,&rhog,P,&err);
90        double dT = (Tc-Tt)/NSAT;
91        i = (int)round(((t - Tt)/(Tc - Tt)*(NSAT-1)));
92        double delt = t - ( Tt + i*dT);
93    
94        double rhofT,rhogT;
95        rhofT =   P->table->satFRho[i] + delt*P->table->satFdRhodt[i];// + 0.5*delt*delt*P->table->satFd2RhodT2[i];
96        rhogT =   P->table->satGRho[i] + delt*P->table->satGdRhodt[i];// + 0.5*delt*delt*P->table->satGd2RhodT2[i];
97    
98        MSG("SAT TEST ::: %f   %f",100*(rhof-rhofT)/rhof,100*(rhog-rhogT)/rhog);
99    
100    /*****************************************Single Phase Table Testing*****************************************/
101      #define NPOINTS 100000      #define NPOINTS 100000
102      double temp_s = 400;      double temp_s = 650;
103      double temp_f = 1200;      double temp_f = 1650;
104      int nT = NPOINTS;      int nT = NPOINTS;
105      double dT = (temp_f-temp_s)/nT;      dT = (temp_f-temp_s)/nT;
   
106    
107    
108      double pressH[NPOINTS],enthalpyH[NPOINTS];      double pressH[NPOINTS],enthalpyH[NPOINTS];
# Line 72  int main(void){ Line 117  int main(void){
117          err = FPROPS_NO_ERROR;          err = FPROPS_NO_ERROR;
118    
119          double T = temp_s + i*dT;          double T = temp_s + i*dT;
120            //FluidState S = fprops_set_Trho(T,rho,P,&err);
121            //pressH[i] = fprops_p(S,&err) ;
122          pressH[i] = P->p_fn(T, rho, P->data,&err) ;          pressH[i] = P->p_fn(T, rho, P->data,&err) ;
123          enthalpyH[i] = P->h_fn(T, rho, P->data,&err) ;          enthalpyH[i] = P->h_fn(T, rho, P->data,&err) ;
124          entH[i] = P->s_fn(T, rho, P->data,&err) ;          entH[i] = P->s_fn(T, rho, P->data,&err) ;
# Line 100  int main(void){ Line 146  int main(void){
146      double msecT = (double)(end - start) / (CLOCKS_PER_SEC/1000);      double msecT = (double)(end - start) / (CLOCKS_PER_SEC/1000);
147    
148    
149      MSG("Percentage Errors");    //  MSG("Percentage Errors");
150      MSG("Temp     \tPressure \tEnthalpy ");    //  MSG("Temp     \tPressure \tEnthalpy ");
151      double pererrp,pererrh;      double pererrp,pererrh;
152      for(i=0; i<nT; ++i){      for(i=0; i<nT; ++i){
153          double T = temp_s + i*dT;          //double T = temp_s + i*dT;
154          pererrp = 100*((pressT[i]-pressH[i])/pressH[i]);          pererrp = 100*((pressT[i]-pressH[i])/pressH[i]);
155          pererrh = 100*((enthalpyT[i]-enthalpyH[i])/enthalpyH[i]);          pererrh = 100*((enthalpyT[i]-enthalpyH[i])/enthalpyH[i]);
156    
# Line 112  int main(void){ Line 158  int main(void){
158      }      }
159    
160    
161      MSG("Percentage Errors");     // MSG("Percentage Errors");
162      MSG("Temp     \t\tEntropy  \t\tU        \t\tG        ");    //  MSG("Temp     \t\tEntropy  \t\tU        \t\tG        ");
163      double pererrs,pererru,pererrg;      double pererrs,pererru,pererrg;
164    
165      for(i=0; i<nT; ++i){      for(i=0; i<nT; ++i){
166          double T = temp_s + i*dT;        //  double T = temp_s + i*dT;
167          pererrs = 100*((entT[i]-entH[i])/entH[i]);          pererrs = 100*((entT[i]-entH[i])/entH[i]);
168          pererru = 100*((intuT[i]-intuH[i])/intuH[i]);          pererru = 100*((intuT[i]-intuH[i])/intuH[i]);
169          pererrg = 100*((gibbsgT[i]-gibbsgH[i])/gibbsgH[i]);          pererrg = 100*((gibbsgT[i]-gibbsgH[i])/gibbsgH[i]);
# Line 141  int main(void){ Line 187  int main(void){
187          ++i;          ++i;
188      }      }
189    
190      MSG("AVERAGE percentage errors for 5 variables p h s u and g");      MSG("AVERAGE percentage errors for 5 variables p h s u and g respectively -->");
191      MSG("%3.6f\t%3.6f\t%3.6f\t%3.6f\t%3.6f\n",av[0]/nT,av[1]/nT,av[2]/nT,av[3]/nT,av[4]/nT);      MSG("%3.6f\t%3.6f\t%3.6f\t%3.6f\t%3.6f",av[0]/nT,av[1]/nT,av[2]/nT,av[3]/nT,av[4]/nT);
192    
193      MSG("Helmholtz did %d calculations in %e seconds", nT*5,msecH/1000);      MSG("Helmholtz did %d calculations in %e seconds", nT*5,msecH/1000);
194      MSG("TTSE did %d calculations in %e seconds", nT*5,msecT/1000);      MSG("TTSE did %d calculations in %e seconds", nT*5,msecT/1000);

Legend:
Removed from v.3023  
changed lines
  Added in v.3024

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