/[ascend]/trunk/models/johnpye/fprops/solve_ph.c
ViewVC logotype

Diff of /trunk/models/johnpye/fprops/solve_ph.c

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

revision 2265 by jpye, Thu Aug 5 08:34:43 2010 UTC revision 2266 by jpye, Thu Aug 5 13:53:01 2010 UTC
# Line 23  Foundation, Inc., 59 Temple Place - Suit Line 23  Foundation, Inc., 59 Temple Place - Suit
23    
24  #include <stdio.h>  #include <stdio.h>
25  #include <math.h>  #include <math.h>
 #include <assert.h>  
26    
27    #define FPE_DEBUG
28    #ifdef FPE_DEBUG
29    # include <assert.h>
30    #else
31    # define assert(ARGS...)
32    #endif
33    
34    #include <setjmp.h>
35    #include <signal.h>
36  #define _GNU_SOURCE  #define _GNU_SOURCE
37  #include <fenv.h>  #include <fenv.h>
38    
39  #define SQ(X) ((X)*(X))  #define SQ(X) ((X)*(X))
40    
41    #ifdef SOLVE_PH_DEBUG
42    # define MSG(STR,...) fprintf(stderr,"%s:%d: " STR "\n", __func__, __LINE__ ,##__VA_ARGS__)
43    #else
44    # define MSG(ARGS...)
45    #endif
46    
47    #define ERRMSG(STR,...) fprintf(stderr,"%s:%d: ERROR: " STR "\n", __func__, __LINE__ ,##__VA_ARGS__)
48    
49    
50  int fprops_region_ph(double p, double h, const HelmholtzData *D){  int fprops_region_ph(double p, double h, const HelmholtzData *D){
51    
52      double Tsat, rhof, rhog;      double Tsat, rhof, rhog;
# Line 47  int fprops_region_ph(double p, double h, Line 64  int fprops_region_ph(double p, double h,
64      return FPROPS_SAT;      return FPROPS_SAT;
65  }  }
66    
67    typedef void SignalHandler(int);
68    
69    jmp_buf mark;
70    void fprops_fpe(int sig){
71        feclearexcept(FE_DIVBYZERO|FE_INVALID|FE_OVERFLOW);
72        longjmp(mark, -1);
73    }
74    
75  int fprops_solve_ph(double p, double h, double *T, double *rho, int use_guess, const HelmholtzData *D){  int fprops_solve_ph(double p, double h, double *T, double *rho, int use_guess, const HelmholtzData *D){
76      double Tsat, rhof, rhog, hf, hg;      double Tsat, rhof, rhog, hf, hg;
77      double T1, rho1;      double T1, rho1;
78      int subcrit = 0;      int subcrit_pressure = 0;
     feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);  
79    
80      if(p < D->p_c){      //feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
81          int res = fprops_sat_p(p, &Tsat, &rhof, &rhog, D);      SignalHandler *old = signal(SIGFPE,&fprops_fpe);
82          if(res){      int jmpret = setjmp(mark);
83              fprintf(stderr,"Unable to solve saturation state\n");      if(jmpret==0){
84              return res;          
85          }          if(p < D->p_c){
86          hf = helmholtz_h(Tsat, rhof, D);              int res = fprops_sat_p(p, &Tsat, &rhof, &rhog, D);
87          hg = helmholtz_h(Tsat, rhog, D);              if(res){
88          fprintf(stderr,"hf = %f kJ/kg, hg = %f\n",hf/1e3,hg/1e3);                  ERRMSG("Unable to solve saturation state");
89                    return res;
90          if(h > hf && h < hg){              }
91              fprintf(stderr,"SATURATION REGION\n");              hf = helmholtz_h(Tsat, rhof, D);
92              /* saturation region... easy */              hg = helmholtz_h(Tsat, rhog, D);
93              double x = (h - hf)/(hg - hf);              MSG("hf = %f kJ/kg, hg = %f",hf/1e3,hg/1e3);
94              *rho = 1./(x/rhog + (1.-x)/rhof);  
95              *T = Tsat;              if(h > hf && h < hg){
96              return 0;                  MSG("SATURATION REGION");
97          }                  /* saturation region... easy */
98                    double x = (h - hf)/(hg - hf);
99                    *rho = 1./(x/rhog + (1.-x)/rhof);
100                    *T = Tsat;
101                    return 0;
102                }
103    
104          subcrit = 1;              subcrit_pressure = 1;
105          if(!use_guess){              if(!use_guess){
106              *T = 1.1 * Tsat;                  *T = 1.1 * Tsat;
107              if(h <= hf){                  if(h <= hf){
108                  *rho = rhof;                      *rho = rhof;
109                  fprintf(stderr,"LIQUID GUESS: T = %f, rho = %f\n",*T, *rho);                      MSG("LIQUID GUESS: T = %f, rho = %f",*T, *rho);
110              }else{                  }else{
111                  *rho = rhog * 0.5;                      *rho = rhog * 0.5;
112                  fprintf(stderr,"GAS GUESS: T = %f, rho = %f\n",*T, *rho);                      MSG("GAS GUESS: T = %f, rho = %f",*T, *rho);
113                    }
114                }
115            }else{
116                if(!use_guess){
117                    *T = D->T_c * 1.03;
118                    *rho = D->rho_c * 0.99;
119              }              }
120          }          }
     }else{  
         if(!use_guess){  
             *T = D->T_c * 1.03;  
             *rho = D->rho_c;  
         }  
     }  
121    
122      fprintf(stderr,"STARTING NON-SAT ITERATION\n");          MSG("STARTING NON-SAT ITERATION");
123      //*rho = 976.82687191126922;          //*rho = 976.82687191126922;
124      //*T = 344.80371310850518;          //*T = 344.80371310850518;
125    
126      T1 = *T;          T1 = *T;
127      rho1 = *rho;          rho1 = *rho;
128      assert(!__isnan(T1));          assert(!__isnan(T1));
129      assert(!__isnan(rho1));          assert(!__isnan(rho1));
130    
131      /* try our own home-baked newton iteration */          /* try our own home-baked newton iteration */
132      int i = 0;          int i = 0;
133      double delta_T = 0;          double delta_T = 0;
134      double delta_rho = 0;          double delta_rho = 0;
135      fprintf(stderr,"STARTING ITERATION\n");          MSG("STARTING ITERATION");
136      while(i++ < 60){          while(i++ < 60){
137          double p1 = helmholtz_p_raw(T1,rho1,D);              double p1 = helmholtz_p_raw(T1,rho1,D);
138          assert(!__isnan(p1));              assert(!__isnan(p1));
139          double h1 = helmholtz_h_raw(T1,rho1,D);              double h1 = helmholtz_h_raw(T1,rho1,D);
140          assert(!__isnan(h1));              assert(!__isnan(h1));
141    
142          fprintf(stderr,"  T = %f, rho = %f\tp = %f bar, h = %f kJ/kg\n", T1, rho1, p1/1e5, h1/1e3);              MSG("  T = %f, rho = %f\tp = %f bar, h = %f kJ/kg", T1, rho1, p1/1e5, h1/1e3);
143          fprintf(stderr,"      p error = %f bar\n",(p1 - p)/1e5);              MSG("      p error = %f bar",(p1 - p)/1e5);
144          fprintf(stderr,"      h error = %f kJ/kg\n",(h1 - h)/1e3);              MSG("      h error = %f kJ/kg",(h1 - h)/1e3);
145    
146          if(p1 < 0){              if(p1 < 0){
147              T1 -= (delta_T *= 0.5);                  T1 -= (delta_T *= 0.5);
148              rho1 -= (delta_rho *= 0.5);                  rho1 -= (delta_rho *= 0.5);
149              continue;                  continue;
150          }              }
151    
152          if(fabs(p1 - p) < 1e-6 && fabs(h1 - h) < 1e-8){              if(fabs(p1 - p) < 1e-6 && fabs(h1 - h) < 1e-8){
153              fprintf(stderr,"Converged to T = %f, rho = %f, in homebaked Newton solver", T1, rho1);                  MSG("Converged to T = %f, rho = %f, in homebaked Newton solver", T1, rho1);
154              *T = T1;                  *T = T1;
155              *rho = rho1;                  *rho = rho1;
156              return 0;                  return 0;
157          }              }
158          /* calculate step */              /* calculate step */
159          double f = log(p1) - log(p);              double f = log(p1) - log(p);
160          double g = h1 - h;              double g = h1 - h;
161          assert(!__isnan(f));              assert(!__isnan(f));
162          assert(!__isnan(g));              assert(!__isnan(g));
163          assert(rho1 != 0);              assert(rho1 != 0);
164          assert(p1 != 0);              assert(p1 != 0);
165          assert(!__isnan(fprops_non_dZdT_v('p', T1,rho1, D)));              assert(!__isnan(fprops_non_dZdT_v('p', T1,rho1, D)));
166          double f_T = 1./p1 * fprops_non_dZdT_v('p', T1,rho1, D);              double f_T = 1./p1 * fprops_non_dZdT_v('p', T1,rho1, D);
167          double f_rho = -1./p1/SQ(rho1) * fprops_non_dZdv_T('p', T1, rho1, D);              double f_rho = -1./p1/SQ(rho1) * fprops_non_dZdv_T('p', T1, rho1, D);
168          double g_T = fprops_non_dZdT_v('h', T1,rho1, D);              double g_T = fprops_non_dZdT_v('h', T1,rho1, D);
169          double g_rho = -1./SQ(rho1) * fprops_non_dZdv_T('h', T1, rho1, D);              double g_rho = -1./SQ(rho1) * fprops_non_dZdv_T('h', T1, rho1, D);
170          assert(!__isnan(f_T));              assert(!__isnan(f_T));
171          assert(!__isnan(f_rho));  
172          assert(!__isnan(g_T));              if(__isnan(f_rho)){
173          assert(!__isnan(g_rho));                  ERRMSG("     rho1 = %f, T1 = %f",rho1, T1);
174                }
175          double det = g_rho * f_T - f_rho * g_T;              assert(!__isnan(f_rho));
176          assert(det!=0);  
177          assert(!__isnan(det));              assert(!__isnan(g_T));
178          fprintf(stderr,"      ∂f/∂T = %e\t\t∂f/∂rho = %e\n",f_T, f_rho);              assert(!__isnan(g_rho));
179          fprintf(stderr,"      ∂g/∂T = %e\t\t∂g/∂rho = %e\n",g_T, g_rho);  
180                double det = g_rho * f_T - f_rho * g_T;
181                assert(det!=0);
182                assert(!__isnan(det));
183                MSG("      ∂f/∂T = %e\t\t∂f/∂rho = %e",f_T, f_rho);
184                MSG("      ∂g/∂T = %e\t\t∂g/∂rho = %e",g_T, g_rho);
185            
186          delta_T = -1./det * (g_rho * f - f_rho * g);              delta_T = -1./det * (g_rho * f - f_rho * g);
187          delta_rho = -1./det * (f_T * g - g_T * f);              delta_rho = -1./det * (f_T * g - g_T * f);
188          assert(!__isnan(delta_T));              assert(!__isnan(delta_T));
189          assert(!__isnan(delta_rho));              assert(!__isnan(delta_rho));
190          fprintf(stderr,"          ΔT   = %f\n", delta_T);              MSG("          ΔT   = %f", delta_T);
191          fprintf(stderr,"          Δrho = %f\n", delta_rho);              MSG("          Δrho = %f", delta_rho);
192    
193          if(subcrit){              if(subcrit_pressure){
194              if(h > hg){                  if(h > hg){
195                  /* vapour */                      /* vapour */
196                  int i = 0;                      int i = 0;
197                  while(rho1 + delta_rho > rhog && i++ < 20){                      while(rho1 + delta_rho > rhog && i++ < 20){
198                      delta_rho *= 0.5;                          delta_rho *= 0.5;
199                        }
200                    }else{
201                        /* liquid */
202                        while(rho1 + delta_rho < rhof && i++ < 20){
203                            delta_rho *= 0.5;
204                        }
205                        if(T1 + delta_T < D->T_t) delta_T = 0.5 * (T1 + D->T_t);
206                        if(T1 < D->T_t) delta_T = +1;
207                  }                  }
208              }else{              }else{
209                  /* liquid */  #if 0
210                  while(rho1 + delta_rho < rhof && i++ < 20){                  /* supercritical... stay above critical temperature of density > rho_crit */
211                      delta_rho *= 0.5;                  if(rho1 + delta_rho > D->rho_c && T1 + delta_T < D->T_c){
212                        delta_T = 0.5 *(T1 + D->T_c);
213                  }                  }
214                  if(T1 + delta_T < D->T_t) delta_T = 0.5 * (T1 + D->T_t);  #endif
                 if(T1 < D->T_t) delta_T = +1;  
215              }              }
216          }else{              /* don't go too dense */
217              /* supercritical... stay above critical temperature of density > rho_crit */              if(rho1 + delta_rho > 2000) delta_rho = 2000;
             if(rho1 + delta_rho > D->rho_c && T1 + delta_T < D->T_c){  
                 delta_T = 0.5 *(T1 + D->T_c);  
             }  
         }  
         /* don't go too dense */  
         if(rho1 + delta_rho > 2000) delta_rho = 2000;  
218                    
219          /* don't go too hot */              /* don't go too hot */
220          if(T1 + delta_T > 5000) delta_T = 5000 - T1;              if(T1 + delta_T > 5000) delta_T = 5000 - T1;
221    
222          /* avoid huge step */              /* avoid huge step */
223          while(fabs(delta_T / T1) > 0.6){              while(fabs(delta_T / T1) > 0.6){
224              delta_T *= 0.5;                  delta_T *= 0.5;
225          }              }
226          while(fabs(delta_rho / rho1) > 0.2){              while(fabs(delta_rho / rho1) > 0.2){
227              delta_rho *= 0.5;                  delta_rho *= 0.5;
228          }              }
229    
230          T1 = T1 + delta_T;              T1 = T1 + delta_T;
231          rho1 = rho1 + delta_rho;              rho1 = rho1 + delta_rho;
232            }
233        }else{
234            /* an FPE occurred */
235            MSG("An FPE occurred");
236      }      }
237    
238      *T = T1;      *T = T1;
239      *rho = rho1;      *rho = rho1;
240        ERRMSG("Iteration failed");
241      return 999;      return 999;
242    
243      int res = fprops_nonsolver('p','h',p,h,T,rho,D);      int res = fprops_nonsolver('p','h',p,h,T,rho,D);
244      fprintf(stderr,"%s: Iteration failed in nonsolver\n",__func__);      ERRMSG("Iteration failed in nonsolver");
245      return res;      return res;
246  }  }
247    

Legend:
Removed from v.2265  
changed lines
  Added in v.2266

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